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Zusammenfassung 


Aufgrund ihrer hohen spezifischen Festigkeit und Steifigkeit bieten 
faserverstärkte Polymere (FRP) ein hohes Leichtbaupotenzial. Die 
gute Formbarkeit diskontinuierlicher (DiCo) FRP erlaubt eine hohe 
Designfreiheit bei vergleichsweise geringen Kosten. Kontinuierliche 
(Co) FRP zeigen außergewöhnliche richtungsabhängige mechanische 
Eigenschaften. Die hier vorgestellte Forschung wurde im Rahmen 
des internationalen Graduiertenkollegs “Integrated engineering of 
continuous-discontinuous long fiber reinforced polymer structures” 
(IRTG 2078) durchgeführt. Hierbei soll DiCoFRP lokal mit CoFRP 
verstärkt werden, so dass der neue Hybridwerkstoff (CoDiCoFRP) die 
individuellen Vorteile optimal kombiniert. 

Betrachtet wird ein Sheet Molding Compound (SMC) Komposit auf 
Basis eines ungesättigten Polyester-Polyurethan-Hybrid (UPPH) Harzes. 
Diese Duromerklasse ermöglicht eine stabile Zwischenkonfiguration 
des Komposits und erlaubt einen integrierten Herstellungsprozess 
von CoDiCoFRP Bauteilen. Der Prozess führt zu einer inhomogenen 
Mikrostruktur mit anisotroper Orientierungsverteilung. Die nachfol- 
gende Arbeit beschäftigt sich mit der mikromechanischen Modellierung 
der Schädigungsentwicklung im Komposit unter Berücksichtigung der 
charakteristischen Faserbündelmikrostruktur. 

Ein auf einer zufälligen sequentiellen Addition basierender Algorithmus 
ermöglicht die schnelle Generierung von SMC Mikrostrukturen. Die 
Verwendung einer exakten Schließung in zwei Dimensionen in Kombina- 
tion mit einer quasi-zufälligen Orientierungsrasterung, führt zu einer be- 
merkenswerten Genauigkeit und ermöglicht die Erzeugung von äußerst 
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realitätsnahen Mikrostrukturen. Dies erlaubt eine Sensitivitatsanalyse 
der effektiven elastischen Eigenschaften des Komposites mit Hinblick 
auf mikrostrukturelle Parameter. Zudem werden die Ergebnisse mit 
direkten numerischen Simulationen auf digitalen Volumenbildern und 
mit “Mean-Field” Methoden verglichen. 

Im Weiteren wird ein Modell zur Berechnung anisotroper Schädi- 
gungsentwicklung im Kontext verallgemeinerter Standardmaterialien 
eingeführt. Unter Ausnutzung der Konvexität der elastischen En- 
ergiedichte in Bezug auf Dehnung und Nachgiebigkeit ist eine thermo- 
dynamisch konsistente Formulierung möglich, die zudem das Schädi- 
gungskriterium nach Wulfinghoff erfüllt. Ein modularer Aufbau unter 
Verwendung von Extraktionstensoren und Schädigungsfunktionen, in 
Analogie zur Elastoplastizität, erlaubt ein breites Anwendungsspektrum 
und die Beschreibung komplexer Schädigungsentwicklungen. Eine 
effiziente Integration der Biot-Gleichung wird diskutiert und die 
Vielseitigkeit des Modells wird an verschiedenen Beispielen auf der 
Mikro- und Makroskala demonstriert. 


Untersuchungen zur Verteilung der Bruchfestigkeit und Steifigkeitsre- 
duzierung von SMC Proben werden ausgewertet. Um die entsprechende 
Schädigungsentwicklung zu erfassen, werden Extraktionstensoren 
eingeführt, die durch die Laminattheorie von Puck motiviert sind. Diese 
bilden verschiedene Schädigungsmechanismen ab und berücksichten 
sowohl gemittelte als auch maximale Spannungszustände. Um die 
Schädigungsparameter zu identifizieren, wird eine Bayes’sche Opti- 
mierung mit Gauß’schen Prozessen verwendet. Abschließend wird 
ein ganzheitlicher Multiskalen-Ansatz zur Identifikation anisotroper 
Versagenskriterien basierend auf Vollfeldsimulationen der mikroskali- 
gen Schädigungsentwicklung präsentiert. Darauf basierend werden 
Versagensflächen im Spannungsraum und durch Steifigkeitsreduktion 
induzierte Versagensflächen vorgestellt, um sowohl die Perspektive 
einer Strukturanalyse als auch die eines Designprozesses abzudecken. 
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Fiber reinforced polymers (FRP) offer lightweight potential due to their 
high specific strength and stiffness. The good formability of discountin- 
uous (DiCo) FRP enables a high freedom in design at comparatively 
low costs in mass production. Continuous (Co) FRP show extraordi- 
nary directional-dependent mechanical properties. The here presented 
research was conducted within the International Research Training 
Group “Integrated engineering of continuous-discontinuous long fiber 
reinforced polymer structures” (IRTG 2078). The development and 
characterization of a new hybrid material combining DiCoFRP with local 
reinforcing CoFRP (CoDiCoFRP) to make the most of their individual 
advantages is the main goal of the IRTG 2078. 

In this thesis, we are concerned with an E-glass fiber reinforced sheet 
molding compound (SMC) composite based on an unsaturated polyester 
polyurethane hybrid (UPPH) resin. This thermoset class enables a stable 
bi-stage condition of the composite and hence allows for a co-molding 
manufacturing process of CoDiCoFRP components. The manufacturing 
process leads to an inhomogeneous microstructure of the SMC com- 
posite with an anisotropic orientation distribution. Our work aims at 
micromechanical modeling of the damage evolution in the composite 
while accounting for the characteristic fiber bundle microstructure. 

We introduce an algorithm based on random sequential addition that 
allows for a fast generation of SMC composite microstructures. Using 
an exact closure approximation in two dimensions in combination with 
a quasi-random orientation sampling derived from a Sobol’ sequence 
leads to a remarkable accuracy and enables the generation of high fidelity 
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microstructures. We comprehensively investigate the sensitivity of the 
effective elastic properties of our SMC composite w. r. t. microstructural 
parameters such as phase properties or volume fraction. We compare the 
results using generated microstructures to direct numerical simulations 
on large scale digital volume images and mean-field estimates. 

We introduce a framework for anisotropic damage evolution in the 
context of generalized standard materials. Exploiting the joint convexity 
of the elastic energy density in terms of strain and compliance, we attain 
the model to be thermodynamically consistent and to fulfill Wulfin- 
ghoff’s damage growth criterion, yielding mesh-independent results. A 
modular formulation of our model using our concept of stress-extraction 
tensors and damage-hardening functions, in analogy to elastoplasticity, 
allows for a wide range of applications and the description of complex 
damage-degradation behavior. We discuss how to efficiently integrate 
Biot’s equation implicitly and demonstrate the versatility of our model 
for a variety of examples on the micro- and macroscale. 

We evaluate experimental investigations regarding distributions of ul- 
timate strength and stiffness reduction of SMC composite specimens. 
To capture the corresponding damage evolution, we introduce stress- 
extraction tensors motivated by Puck’s laminate theory accouting for 
different damage mechanisms in both an averaged and a maximum 
stress setting. To identify all damage parameters associated to our SMC 
composite, we utilize Bayesian optimization with Gaussian processes. 
With these tool at hand, we present a holistic multiscale approach for 
constructing anisotropic failure criteria based on full-field simulations 
of microscale damage evolution. We propose failure surfaces in stress 
space and stiffness-reduction triggered failure surfaces to cover both a 
structural analysis and a design process perspective. 
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Chapter 1 


Introduction’ 


1.1 Motivation 


Weight savings. One key element to combat climate change and its 
impacts is a drastic lowering of carbon dioxide emissions and a long- 
term systemic shift (United Nations Environment Programme, 2015; 
2019). Regarding the automotive sector, the fuel consumptions of cars 
and trucks can be reduced by about 0.40 1 to 0.49 1 per 100 km for each 
100 kg of weight savings. Consequently, 1.20 US$ to 13.70 US$ are being 
paid per kilogram of weight savings in that industry (Bandivadekar 
et al., 2008). 

Fiber reinforced polymers. Lightweight design plays a significant role 
in the reduction of fuel consumption and eventually the carbon dioxide 
emissions. One approach is the substitution of (non-structural) metallic 
parts with lighter fiber reinforced polymer (FRP) components. Generally, 
these composites consist of a matrix resin system, reinforced by fibers 
and conceivably some additives and fillers. Current fiber materials 
include glass, wood, carbon and cotton, see Fig. 1.1b. Typically, matrix 
resins are categorized in thermosets (TS) and thermoplastics (TP). Im- 


1 This chapter is based on excerpts of the publications “Computational homogenization 
of SMC composites based on high fidelity representative unit cells” (Görthofer et al., 
2020), “A convex anisotropic damage model based on the compliance tensor” (Gérthofer 
et al., 2022b) and “A computational multiscale model for anisotropic failure of SMC 
composites” (Görthofer et al., 2022a) 
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proved manufacturing technologies allow for cost-effective productions 
and consequently applications that are not restricted to the automotive 
or aerospace sector, but also include construction, sports or leisure, see 
Fig. 1.1a. 


Sr 


» 


E Transportation 35% E Glass 85% 
E Construction and infrastructure 34% E Wood 7% 
E Sports and leisure 15 % Cotton 4% 
E Electrical and electronic equipment 15% Carbon 2% 
E Others 1% E Others 2% 
(a) Market segments for glass fiber rein- (b) Different fiber types (Reddy, 2015) 


forced composites (Witten et al., 2017) 


Figure 1.1: Distributions of composites and fiber types in Europe (Kehrer, 2019) 


Types of FRP. Depending on the composition, FRP exhibit a high mass 
specific strength and stiffness. Generally, FRP are sub-divided into 
different classes, i.e., discontinuous FRP (DiCoFRP) where (short or 
long) fibers are dispersed more or less randomly and continuous FRP 
(CoFRP) with unidirectionally aligned fibers, see Fig. 1.2. The latter 
enable a loading case specific manufacturing of high performance com- 
ponents due to designated arrangements of the continuous fibers. A 
fiber volume fraction of up to 60 % further supports high strength and 
stiffness properties. The known fiber directions facilitate the modeling 
and dimensioning process of CoFRP components. Drawbacks of this 
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material class include long cycle times, high scrap rates and a limited 
formability. Complementary, DiCoFRP generally show lower stiffness 
and strength properties compared to CoFRP, but offer a high freedom 
in design at short cycle times. The manufacturing processes, including 
compression molding and injection molding, allow for an economic high- 
volume production of complex components, counting, e. g., beads or ribs. 
DiCoFRP are currently applied as non-structural components (Ernst 
et al., 2006), and thus make a first contribution towards lightweight 
applications as replacements for certain metallic parts. As the complex 
mold-filling process leads to a process-dependent and heterogeneous 
fiber orientation distribution, the development of efficient modeling 
approaches and dimensioning tools is a challenging task. 


CoFRP CoDiCoFRP DiCoFRP 


Continuous Fiber Discontinuous Fiber Reinforced Polymer Discontinuous Fiber 
Reinforced Polymer with Continuous Fiber Reinforcement Reinforced Polymer 


Figure 1.2: Classification of different fiber reinforced polymers 


Hybridization. A combination of CoFRP and DiCoFRP to a new hybrid 
class of continuous-discontinuous FRP (CoDiCoFRP), see Fig. 1.2, aims 
at merging the individual advantages. Specific reinforcements of CoFRP 
within the DiCoFRP part are used to sustain application-based loadings 
where necessary. Hence, CoDiCoFRP offer a high potential for the 
production of composite parts to be used as structural components with 
a high specific strength and stiffness at comparably low cycle times and 
reduced costs. A CoDiCoFRP demonstrator part is shown in Fig. 1.3, 
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where the black areas are carbon fiber CoFRP reinforcements within the 
lighter glass fiber DiCoFRP. 


Figure 1.3: IRTG 2078 demonstrator part made of glass fiber DiCoFRP with local carbon 
fiber CoFRP reinforcements (manufactured at Fraunhofer Institute of Chemical Technology 
(ICT) Pfinztal) 


IRTG 2078. In order to unlock the full potential of 
these hybrid materials, a profound understand- 
ing, covering technological manufacturing and 
experimental characterization aspects, as well as 
component design and simulation prospects, is 


ERP 


of importance. Therefore, the complex multiscale 
material is analyzed in an integrated engineering Figure 1.4: Logo of 
approach within the International Research Train- IRTG 2078 

ing Group "Integrated engineering of continuous- 

discontinuous long fiber reinforced polymer struc- 

tures" (IRTG 2078). The international and interdisciplinary collaboration 
links necessary expertise to propel research on CoDiCoFRP aiming 
at a full understanding of the material’s behavior and development 


of associated simulation methods. To fully cover both the physical 
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and virtual process chain for CoDiCoFRP components, IRTG 2078 is 
sub-divided into four research areas, i. e., technology, characterization, 
design and simulation, see Fig. 1.5. 


Technology Characterization 


Preforming Molding 
A 
v 


Concept Optimize Prediction Microstructure Damage 


ro 


Simulation 


Figure 1.5: Schematic of physical and virtual process for CoDiCoFRP (Görthofer et al., 
2019a) including an interface CT scan performed at Fraunhofer Institute for Mechanics of 
Materials (IWM) Freiburg (Schober et al., 2017) 


Physical process chain. The physical process chain is established within 
the research areas technology and characterization. Key aspects are 
the application-based preforming of CoFRP tapes (Kupzik et al., 2020), 
the simultaneous co-molding of DiCoFRP and CoFRP while ensuring 
proper interface bonding (Biicheler, 2018) and the subsequent machin- 
ing processes (Langer et al., 2020). A continuous surveillance of the 
individual production steps is important to ensure a consistent high 
quality of the manufactured parts (Bretz et al., 2021). Accompanying 
the manufacturing process, experimental investigations are conducted 
ranging from microscale analysis to macroscopic structural component 
testing. Interface characterization provides information about fiber 
and matrix bonding (Rohrmiiller et al., 2020). Local variations of fiber 
orientations are identified via micro-computed tomography analysis 
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(Schöttl et al., 2021b). Component scale testings are used to identify 
effective properties of DiCoFRP and CoFRP and the reinforcing effects 
of the hybrid material regarding, inter alia, quasi-static (Trauth et al., 
2017a) and fatigue loading scenarios (Bartkowiak et al., 2020). 


Virtual process chain. The associated virtual process chain mimics 
the physical process chain providing a constant data and knowledge 
exchange to improve the manufacturing of CoDiCoFRP components. 
Within research area design a consistent system of objectives is estab- 
lished to gather all achieved progress and provide structured concept 
criteria for the working engineer (Richter et al., 2020). CoDiCoFRP 
components are optimized w.r.t. their shape, usage of beads and ribs, 
as well as positioning of CoFRP tapes to maximize application-based 
criteria such as strength or stiffness (Fengler et al., 2019; Revfi et al., 2021). 
Process simulation techniques are developed to correctly predict the flow 
of the DiCoFRP during molding and hence allow for the prediction of 
inhomogeneous orientation states within the manufactured parts (Meyer 
et al., 2020). In order to describe the structural behavior of CoDiCoFRP 
components, different key aspects are analyzed in the research area 
simulation. Curing simulations on the microscale help quantifying 
residual stresses (Pallicity and Böhlke, 2021; Schwab, 2019). Mean- 
field (Schemmann et al., 2018b) and full-field models (Görthofer et al., 
2022b) taking into account the FRP microstructure, phase properties 
and different damage mechanisms allow for an micromechanics-based 
evaluation of the effective stiffness reduction. The junction of these 
simulative approaches enables a component scale computation of the 
structural behavior (Gorthofer et al., 2019b). 


1.2 Objectives and originality of this thesis 


Categorization. In this thesis we present investigations and results 
conducted within the scope of project “S3: micro-mechanical modeling” 
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in research area simulation of the second generation in the IRTG 2078. 
The material system under consideration is sheet molding compound 
(SMC), which is a thermoset (TS) matrix reinforced with glass and 
carbon fibers within the DiCo and Co phases, respectively. To ensure a 
proper co-molding of this CoDiCoFRTS hybrid, an unsaturated polyester 
polyurethane hybrid (UPPH) resin is used (Bücheler et al., 2017). 
Objectives. The main objective of our presented research is the devel- 
opment of a micromechanical model, incorporating damage, that helps 
to propel scale transition and provide micromechanics-based input for 
the macroscopic component scale simulations. Therefore, we mainly 
address three topics in this thesis. 

(I) Ina first step, we establish a microstructure generator that ac- 
counts for the characteristic fiber bundle structure of SMC com- 
posites, see Fig. 1.6, and enables the generation of high fidelity 
representative unit cells. The means to individually adjust any 
volume fraction and orientation state within the unit cells serves 
as basis for a sensitivity analysis of effective elastic properties, 
as well as a utilization of the generated microstructures in the 
further process of material modeling (Görthofer et al., 2020). 
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In a second step, we propose a general damage model framework 
that allows for a description of any anisotropic damage evolution 
and associated stiffness reduction. The model is formulated in 
the setting of generalized standard materials (GSM) of dissipative 
solids, ensuring thermodynamical consistency. A convex mod- 
eling of the free energy density and the force potential enables a 
unique solution for any material that exhibits a damage hardening 
regime prior to total failure. The compliance tensor serves as 
natural and observable internal damage variable (Görthofer et al., 
2022b). 

(III) In the third step, we apply our presented model to SMC com- 
posites in order to predict damage evolution on the microscale. 
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Therefore, we identify corresponding damage cases, inspired 
by Puck’s criteria for laminates. We design specific extraction 
tensors that describe matrix and fiber bundle damage in SMC 
composites, taking into account either an averaged stress state 
or maximum stresses. We use a Bayesian optimization approach 
to identify all associated SMC composite damage parameters. 
Using our generated unit cells, we compute a variety of loading 
scenarios on different microstructures to analyze the damage 
evolution. A comparison to pCT scans and experimental results 
yields micromechanics-based anisotropic failure criteria to be 
used on a component scale level (Görthofer et al., 2022a). 


Originality of this thesis. Our presented research on SMC composites 


and damage modeling exhibits the following novelties: 


SMC composite characteristics: We explicitly consider the inherent 
three-scale structure of SMC composites taking into account the char- 
acteristic fiber bundles. Furthermore, we account for the macroscopic 
strain hardening behavior due to accumulation of damage on the 
microscale that eventually leads to abrupt failure. 

Microstructure generation: Using a quasi-random orientation sam- 
pling based on a Sobol’ sequence in combination with the exact closure 
approximation allows for an accurate microstructure reconstruction 
yielding high fidelity representative SMC composite unit cells. We 
present an algorithm that enables a fast and robust generation of SMC 
composite microstructures based on random sequential addition. 
Sensitivity analysis: Utilizing a large variety of generated unit cells, 
we provide a comprehensive sensitivity analysis of effective stiffness 
properties w.r.t. changing phase properties, volume fraction and 
orientation. 

Anisotropic and convex damage framework: Our anisotropic dam- 
age model is not only derived in the GSM framework, thermody- 
namically consistent and features a convex free energy and a convex 
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dissipation potential, but also fulfills Wulfhinghoff’s damage growth 
criterion. Our presented framework is local, but nonetheless well- 
posed and hence leads to a unique, mesh-independent solution. We 
show and make use of the fact that the elastic energy density of linear 
elasticity is jointly convex in the strain and the compliance tensor. 
The compliance tensor as a primary damage variable serves as an 
experimentally measurable internal state variable. 

e Modularity: A modular framework in combination with extraction 
tensors enables the application of our model to a broad set of materials 
on the micro- and the macroscale. We introduce specific extraction 
tensors motivated by Puck’s laminate theory. Based on an averaging 
process, as well as a pencil glide approach, we define both extraction 
tensors accounting for average stress states and maximum stresses, re- 
spectively. Furthermore, we introduce loading case specific extraction 
tensors accounting, inter alia, for distortion and dilatation. 

e Micromechanis-based failure criteria: We present a holistic multi- 
scale approach for constructing anisotropic criteria describing the 
macroscopic failure based on microscale full-field damage evolution. 
Failure surfaces in stress space and stiffness-reduction triggered fail- 
ure cover both a structural analysis and a design process perspective. 


1.3 State of the art 


1.3.1 Research on SMC composites 


Characteristics. For lightweight applications, DiCoFRP combine low 
material and manufacturing costs with a high degree of freedom in 
design (Wilkinson and Ryan, 1998). In particular, sheet molding com- 
pound (SMC) composites are frequently used in applications due to 
their high strength-to-weight ratio (Huang and Zhao, 2012; Asadi et al., 
2017). Standard SMC composite products consist of discontinuous fiber 
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bundles, a mixture of thermosetting resin (typically, polyester, vinylester 
or epoxy), fillers (calcium carbonate, alumina, etc.) and additives (such 
as initiators, inhibitors, thickeners), see Dumont et al. (2007) or Kim 
et al. (2011a). The compression molding manufacturing process allows 
for producing components at comparatively low costs while offering a 
high freedom in design (Wilkinson and Ryan, 1998). Furthermore, the 
manufacturing process of SMC leads to a characteristic microstructure 
where fibers are almost aligned in bundles (Le et al., 2008; Meyer et al., 
2020; Schöttl et al., 2021b), see Fig. 1.6. 


Microscale Mesoscale Macroscale 


Figure 1.6: Three-scale structure of SMC composites, see Görthofer et al. (2019b) 


Hybridization. Combining glass and carbon fiber bundles (Brinson and 
Brinson, 2008; Anagnostou et al., 2018; Trauth, 2020) permits applying 
SMC composites as load-bearing structural components. A modern 
unsaturated poylester polyurethane hybrid (UPPH) resin without fillers 
allows to manufacture SMC composites with discontinuous glass fiber 
bundles and local continuous carbon fiber patch reinforcements within 
a simultaneous co-molding process (Biicheler et al., 2017; Bohlke et al., 
2019). 

Experimental investigations. The increased use of SMC composites 
fosters the development of accurate mechanical models (Görthofer et al., 
2019b; Böhlke et al., 2019). Experimental investigations of the microstruc- 
ture (Motaghi and Hrymak, 2017; Schöttl et al., 2020; 2021b), the effective 
elastic behavior in a quasi-static (Irauth et al., 2017a; Chen et al., 2018b) 
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and dynamical setting (Fitoussi et al., 2005; Jendli et al., 2005; Kehrer 
et al., 2018) were performed. In-situ testing using micro-computed 
tomography (nCT) sheds light on the underlying damage mechanisms 
(Dumont et al., 2007; Arif et al., 2014; Rohrmiiller et al., 2020) and 
microcrack distributions (Schöttl et al., 2020; Schöttl et al., 2021a). 

The influence of temperature, especially with an eye towards the glass 
transistion of the polymer resin, was analyzed by Kehrer et al. (2018) 
for the elastic regime. Meyer, Hohberg and co-workers studied the flow 
and curing behavior during the manufacturing process (Hohberg et al., 
2017; Meyer et al., 2020). The behavior of SMC composites under fatigue 
loading was studied by Ben Cheikh Larbi et al. (2006) and Bartkowiak 
et al. (2019; 2020). Investigations on the structural behavior of SMC 
composites under combined stress states were performed via dedicated 
cruciform specimens (Ogihara and Koyanagi, 2010; Schemmann et al., 
2018a;c). Non-destructive measurement techniques were developed for 
quantifying manufacturing uncertainties and hence the quality of SMC 
composite parts (Schäferling et al., 2018; Bretz et al., 2020; 2021). 
Damage and failure. Further focus of recent research on SMC com- 
posites is devoted to damage and failure. Fitoussi et al. (1996; 1998) 
presented approaches for predicting anisotropic damage evolution based 
on multilocal criteria, also at high strain rate (Fitoussi et al., 2013). Matrix 
degradation and interface decohesion, two major damage mechanisms 
in SMC composites, were analyzed and modeled by Meraghni and 
co-workers (Meraghni and Benzeggagh, 1995; Meraghni et al., 1996) 
and implemented in a micromechanical framework (Meraghni et al., 
2002). A similar approach for predicting the non-linear behavior of 
SMC composites was proposed by Baptiste (2003). Dedicated models 
were developed that account for the multiscale and anisotropic nature 
of SMC composites (Drugan and Willis, 1996a; Guo et al., 1997), taking 
into account reliability (Desrumaux et al., 2000), inclusion distributions 
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(Zheng and Du, 2001; Lee and Simunovic, 2001), humidity (Arif et al., 
2014) or dynamical behavior (Jendli et al., 2009). 

SMC composites typically show an elasto-damageable behavior that is 
driven by anisotropic damage evolution (Fitoussi et al., 2005; Ogihara 
and Koyanagi, 2010) and eventually ends in abrupt failure. Damage 
evolution on the microscale induces an effective stress-strain relationship 
characterized by a hardening, rather than a softening, regime (Trauth 
et al., 2017a; Trauth, 2020). Hence, no localization occurs prior to a 
specific loading point right before total failure of the SMC composite 
(Meraghni and Benzeggagh, 1995; Görthofer et al., 2019b). 


1.3.2 Microstructure generation 


Motivation. To reduce safety factors in dimensioning SMC parts, ac- 
curate material models need to be provided, accounting both for the 
intrinsic complexity of the composite behavior - temperature- and strain- 
rate dependent mechanical behavior involving damage (Schemmann 
et al., 2018b; Shirinbayan et al., 2017) - and the manufacturing process 
induced variation of microstructural parameters (such as fiber volume 
content and fiber orientation). 

To account for the inherent anisotropy of the effective material behavior 
of fiber reinforced composites, both analytical (Mori and Tanaka, 1973; 
Ponte Castafieda and Suquet, 1998; Duschlbauer et al., 2003; Doghri and 
Tinel, 2005; Anagnostou et al., 2018) and computational homogenization 
techniques (Matous et al., 2017) are frequently used to obtain accurate 
effective material models and to identify the physical mechanisms re- 
sponsible for the effective mechanical behavior. 

SMC composite unit cells. Using representative unit cells as the basis 
for computational homogenization has become a standard technology 
for short fiber reinforced composites, cf. the overview articles of Ma- 
tous et al. (2017) and Bargmann et al. (2018). Of key importance here 
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is generating high fidelity volume elements, accounting both for the 
high fiber volume content and the fiber orientation to high precision 
(Schneider, 2017). For long fiber reinforced composites, the fibers can no 
longer be described by straight cylinders, and more elaborate methods 
for generating unit cells need to be employed, see for instance Altendorf 
and Jeulin (2011) or Fliegener et al. (2014). 

Due to differences in the production 
process, the microstructures of long 
fiber reinforced thermoplastics (LFT) 
and SMC composites exhibit different 
characteristics. In contrast to the for- 
mer, distinct fiber bundles are appar- 
ent for SMC composites. In particular, 
unit cell generation techniques for 
LFT cannot be applied to SMC in a 
straightforward fashion. For illus- 


tration, Fig. 1.7 shows a microstruc- 

ture generated by the SAM algorithm Figure 1.7: LFT microstructure gener- 

(Schneider, 2017) for a fiber orienta- ated by the SAM algorithm (Schneider, 
. : 3 i 2017) for 22.5 % fibers with an aspect ra- 

tion that is almost planar isotropic tio of 100 and a second-order fiber orien- 

and fibers with an aspect ratio of 100. tation tensor of diag(0.475, 0.475, 0.05), 

. _ see Görthofer et al. (2020) 

The edge length of the cubical unit 

cell is roughly twice the fibers’ length. 

For a fiber volume fraction of 22.5 %, the cell contains 24682 fibers. For 

this kind of microstructure, computational homogenization is difficult. 

Furthermore it should be noted that the aspect ratio for the fibers used 

in SMC composites is even one order of magnitude higher. Furthermore, 

the induced bundle structure of the LFT microstructure in Fig. 1.7 is 


different from the bundle structure observed in SMC, cf. Fig. 1.6. 


Physical approaches. For this reason, unit cell modeling of SMC com- 
posites is focused on the apparent three-scale structure, as shown in 
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Fig. 1.6, modeling fiber bundles (sometimes also called tows) directly. 
Extending the work of Fliegener et al. (2014) on LFT, Fette et al. (2017) 
used a commercial finite element code to compress previously placed 
fiber bundles accounting for developing contact between the fiber bun- 
dles. A similar approach was presented by Li et al. (2018b). Islam et al. 
(2016) coupled a random sequential addition algorithm to a dynamic 
finite element simulation. Sparse inclusion assemblies are generated in 
pseudo-boxes. A subsequent explicit flow simulation is used to transfer 
the inclusions from the pseudo-boxes into the target box. For all of these 
approaches, constructing these volume elements is computationally 
most demanding, and realizing a prescribed fiber orientation state is 
based on a trial-and-error procedure. Altendorf and Jeulin (2011) use 
a force-biased packing approach. Fibers are represented as chains of 
balls with a controllable level of bending. A deposition algorithm in 
combination with a force-directed approach based on Coulomb’s law 
and spring forces to hinder bundle overlap was presented by Harper et al. 
(2017). Smooth fiber undulations are ensured via spline interpolations. 

Geometric approaches. A periodic packing of general ellipsoids was 
presented by Ghossein and Lévesque (2013), where ellipsoids can move, 
rotate and collide in order to reach a desired volume fraction. Li et al. 
(2018a) studied a Voronoi diagram based algorithm for modeling SMC 
microstructures, where each Voronoi cell corresponds to a fiber bundle 
with a distinct orientation. These Voronoi cells recover the characteristic 
tesselated surfaces of carbon fiber SMC with high fiber volume content. 
As typical for Voronoi-based modeling, it is difficult to realize a given 
statistics of the cell geometry (Quey and Renversade, 2018). Motivated 
by the random sequential adsorption algorithm (Feder, 1980), Chen et al. 
(2018b) developed an algorithm for modeling SMC unit cells based on 
a sequential bundle deposition method. Chen et al. (2018b) model the 
SMC fiber bundles as cuboid cells which are placed sequentially and 
layer by layer, accounting for an overlap between bundles by displacing 
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this overlap to an adjacent layer. The advantage of this algorithm 
is that both the bundle dimensions and the fiber orientation can be 
prescribed independently. Furthermore, Chen et al. (2018b) showed a 
close agreement of their computational results with experimental data 
for carbon fiber SMC. 


1.3.3 Damage modeling 


Progressive degradation. Damage mechanics describes the progres- 
sive degradation of the elastic stiffness of materials upon loading, and 
is typically attributed to growing voids or cracks on a lower length 
scale (Lemaitre, 1996), see Fig. 1.8. There are two predominant ap- 
proaches to continuum damage-mechanics (Krajcinovic, 1984; Lemaitre 
and Chaboche, 1990). The first approach accounts for the origin of 
damage on a lower length scale in terms to micromechanics (Fitoussi 
et al., 1996; Guo et al., 1997), see also Sec. 3 in Krajcinovic (1989) for an 
early account. With qualitative predictions in mind, the second strategy 
is of phenomenological nature. After selecting a suitable damage vari- 
able (or a collection thereof), suitable kinetic laws are postulated taking 
continuum thermodynamics into account, Sec. 4 in Krajcinovic (1989). 


Grains Matrix Fibers 


ON D aX 
Microcracks — Microcracks 


(a) Polycrystalline microstructure (b) Microstructure composed of fiber bundles 
Figure 1.8: Schematics of microstructures with growing microscopic cracks, passing from 


state @ to state @, similar to Fassin et al. (2019). Growing microcracks induce a reduction 
of the effective stiffness. 
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Micromechanics-based approaches. The micromechanics-based ap- 
proach to damage mechanics takes the damage mechanisms on a lower 
scale into account and is still subject of current research, for instance 
concerning mesh-size objective modeling (Liang et al., 2018), a coupling 
to model-order reduction (Bhattacharyya et al., 2020) or accounting for 
micro-computed tomography data (Luo et al., 2020). Micromechanics- 
informed damage models permit taking the stochastics on the microscale 
into account naturally, e. g., for progressive fiber breakage in fiber re- 
inforced composites (Ju and Wu, 2016; Wu and Ju, 2017), interfacial 
transition-zone effects (Chen et al., 2018a), uncertainty in the elastic mod- 
uli of fiber reinforced concrete (Liu et al., 2020), localized microcracks 
(Li et al., 2020) or random loading in fatigue processes (Franko et al., 
2017). Another advantage concerns modeling the unilateral character 
of brittle damage, i.e., a different damaging behavior under tension 
compared to compression (Goidescu et al., 2015; Zhang et al., 2019), and 
accounting for interface debonding (Pupurs and Varna, 2017; Schem- 
mann et al., 2018b; Yang et al., 2020). However, care has to be taken as 
homogenization and localization are incompatible (Gitman et al., 2007), 
in general, i. e., upon localization, the volume elements considered will 
not be representative for the effective mechanical behavior (Hill, 1963; 
Drugan and Willis, 1996b; Kanit et al., 2003). 

Many micromechanics models are based on mean-field methods (Hashin 
and Shtrikman, 1962; Mori and Tanaka, 1973; Willis, 1981), e. g., taking 
into account a numerically computed Eshelby’s tensor (Desrumaux et al., 
2001), an evaluation of inclusion stresses (Duschlbauer et al., 2003), a 
Weibull probability density for the interface strength (Schemmann et al., 
2018b) or a mixture of glass and carbon fibers (Anagnostou et al., 2018). 
Computational homogenization techniques (Schneider, 2021) were used 
for studying the effective elastic behavior of SMC composites (Matouš 
et al., 2017; Görthofer et al., 2020). 
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Phenomenological approaches. As an alternative to micromechanics- 
type strategies, phenomenological approaches to continuum-damage 
mechanics may be pursued. In a first step, a (scalar- or tensor-valued) 
damage variable is selected which describes the reduction of the effective 
cross-section of a typical material sample undergoing material degra- 
dation (Gurson, 1977; Voyiadjis, 2015). Then, suitable kinetic laws are 
postulated on the basis of continuum thermodynamics (Simo and Ju, 
1987; Hansen and Schreyer, 1994). 

Tensor order of damage variable. The tensor order of the damage 
variable naturally distinguishes different phenomenological damage 
models. Even today, the classical scalar isotropic damage variable serves 
as a reliable workhorse with numerous applications including cast steel 
with pores (Yan et al., 2020), concrete (Li and Wu, 2018), rocks (Liu et al., 
2018; Xu et al., 2018), framed structures (Yang et al., 2017), unidirectional 
glass fiber reinforced plastic composite plies (Sharma and Daggumati, 
2020), fibrous composite laminae (Abu-Farsakh and Asfa, 2020), notched 
epoxy resin specimens (Rahimi et al., 2020) and steel-fiber reinforced 
concrete (Moradi et al., 2020). 

Second-order tensors. Damage variables with higher tensor order 
permit modeling an emerging anisotropy of damage. As working with 
second-order tensors comes naturally to disciples of continuum mechan- 
ics, it is not surprising that second-order damage-tensors (Murakami and 
Ohno, 1981) are used frequently in continuum damage-mechanics. Re- 
cent applications include concrete (Desmorat, 2016; Wardeh and Toutanji, 
2017), metal-forming processes (Nasab and Mashayekhi, 2019), rock ma- 
terials (Wang and Xu, 2020), composite fabrics and laminated panels (Wei 
et al., 2020) and composite laminates (Okabe et al., 2018; Onodera and 
Okabe, 2020). Second-order damage-tensors are always orthotropic w. r. t. 
their eigenbasis, limiting their degree of generality. More often than not, 
such a limitation is interpreted as a feature, and specific orthotropic 
damage models are developed, for instance for brittle materials (Kim 


17 


1 Introduction 


et al., 2016), in elastoplastic and finite-strain damage coupling (Ganjiani, 
2018; Reese et al., 2021), or for ceramic-matrix composites (Alabdullah 
and Ghoniem, 2020). 

Fourth-order tensors. As continuum damage-models primarily seek 
to describe a loss of stiffness due to emerging defects in solids, using a 
fourth-order damage-tensor (Chaboche, 1981), the same tensor order as 
the stiffness tensor, appears reasonable. In Sec. 4.3.4, Krajcinovic (1989) 
even notes that "an appropriate description of damage [...] must involve 
at least a fourth-rank tensor." This idea was pursued for the stiffness or 
compliance tensors as the primary damage variable (Dougill, 1976; Ortiz 
and Popov, 1982; Ortiz, 1985), also coupled to plasticity (Simo and Ju, 
1987; Ju, 1989; Yazdani and Schreyer, 1990). We refer to Zhang and Cai 
(2010) for a modern account of anisotropic damage mechanics. However, 
some care has to be taken when working with tensor-valued damage 
variables due to possible inconsistencies arising for complex non-radial 
loading-unloading scenarios, see Simon et al. (2017). 
Tension-compression asymmetry. The unilateral character of pores 
and cracks (see Fig. 1.8) often leads to a tension-compression asymmetry 
of the material behavior upon damage loading, see Chaboche (1993) 
for a discussion. To incorporate the latter effect in phenomenological 
models, one may introduce different damage variables for the tensile 
and the compressive regime (Ramtani et al., 1992; Cicekli et al., 2007). 
For three-dimensional stress states, spectral decompositions of either 
the strain or the stress tensor may form the basis of continuum damage 
models that differentiate between damage evolution due to tension and 
compression (Ladevéze and Lemaitre, 1984; Ortiz, 1985). 

Softening. Whenever damage models exhibit a softening behavior, 
their use in a continuum formulation leads to an ill-posed problem due 
to localization effects (Lemaitre, 1986), which is reflected by strongly 
mesh-dependent results in numerical simulations (De Borst, 1996). Coun- 


termeasures in the framework of local damage models were investi- 


18 


1.3 State of the art 


gated (Tvergaard, 1982; Becker et al., 1988; Beremin et al., 1983). Non- 
local formulations (Belytschko et al., 1986; Bazant, 1991) prevent the 
localization responsible for the ill-posedness, and may be realized by 
an explicit convolution with a tapering function (Pijaudier-Cabot and 
Bazant, 1987), by augmenting the damage evolution equation by an 
elliptic differential operator (Aifantis, 1984) or by employing a gradient- 
enhanced formulation (Briinig and Ricci, 2005; Germain et al., 2007; Abu 
Al-Rub and Voyiadjis, 2009)), which may also be coupled to Hamilton’s 
least-action principle (Junker et al., 2019; 2021). As long as the softening 
is not too pronounced, existence of results for non-local damage models 
(Thomas and Mielke, 2010) may be established. However, except for 
specific models (Susu, 2017; Roubíček, 2009), uniqueness (and, thus, 
well-posedness) cannot be ensured. For a review on ill-posedness due 
to localization problems and appropriate regularization methods, the 
reader is referred to Forest et al. (2004). Also, for a general overview on 
continuum damage-mechanics and further literature, the reader may 
consult the books of Murakami (2012) and of Voyiadjis (2015). 

Hardening. Oftentimes, the ill-posedness of local damage models is 
taken for granted, and appropriate countermeasures are taken. A charm- 
ing strategy takes a conventional local damage model with softening 
(but sufficient growth at infinity), and applies relaxation techniques 
(Balzani and Ortiz, 2012; Schmidt and Balzani, 2016; Schwarz et al., 2021), 
which are typically used for studying solids with emerging microstruc- 
ture. When describing stable damage processes, these countermeasures 
should not be necessary, however. Indeed, for a moderate degree of 
loading, localization is excluded, and manifests only at a specific turning 
point in loading level. For component-scale simulations, this loading 
level is not readily apparent, and depends on the specimen geometry via 
solving the equations of continuum mechanics. To sum up and loosely 
speaking, we know that local damage models are perfectly reasonable 
up to a specific level of loading, but we do not know this level in 
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advance. Thus, interest arose to design damage models which give 
rise to a meaningful response for the entire range of loading, and which 
are intended to be complemented by a classical failure criterion. 


1.4 Outline 


Following the introduction and classification of our research in Chap- 
ter 1, we sub-dived this thesis into two further chapters that present 
some important fundamentals, three chapters based on published jour- 
nal publications discussing the main research focus and a summary 
chapter highlighting the key achievements and possible scopes of further 
research. 

In Chapter 2 we briefly introduce the necessary fundamentals of con- 
tinuum mechanics regarding kinematics and balance equations in the 
context of small deformations. 

Chapter 3 is dedicated to SMC composites, i.e., we present the com- 
pression molding manufacturing process, the resulting characteristic 
microstructure and mathematical descriptions of such. 

In Chapter 4 we present a method that allows the generation of high 
fidelity microstructures of SMC composites via an adapted random 
sequential addition. Using a quasi-random orientation sampling and 
an exact closure, we are able to generate a wide range of unit cells with 
high accuracy. Using full-field homogenization, we conduct a sensitivity 
analysis of effective parameters based on microstructure variations and 
changes of phase properties and compare our results to mean-field 
calculations. 

We introduce a basic framework to model anisotropic damage evolution 
in the context of generalized standard materials in Chapter 5. The 
modular setting allows for a straight-forward adaption of the model to 
different materials using extraction tensors. We discuss different such 


20 


1.5 Notation, frequently used acronyms, symbols and operators 


extraction tensors and subsequently show the capacity and performance 
of our model via different computational examples. 

In Chapter 6 we analyze experimental investigations on SMC composites 
provided by M. Bartkowiak, A. Trauth and L. Schöttl. We utilize our 
damage model and introduce specific extraction tensors motivated by 
Puck’s theory on laminates to capture damage in SMC composites. A 
Bayesian optimization approach is presented to determine all necessary 
damage parameters. Consequently, we present an holistic approach to 
identify anisotropic failure criteria for SMC composites based on the 
computed full-field damage evolution. 

We briefly summarize our work and give concluding remarks and 
conceivable follow-up research in Chapter 7. 


1.5 Notation, frequently used acronyms, sym- 
bols and operators 


We follow a direct tensor notation throughout the text, representing 
vectors and tensors by their components or using matrix representations 
(in an orthonormal basis) only when necessary. Vectors and second- 
order tensors are denoted by lower case and upper case bold letters, 
respectively (e.g., a and A). Fourth-order tensors are denoted by, 


e.g., A, B. Scalars and arrays of quantities are represented by non-bold 
letters (e.g., H, w or z). The transposition of a vector and second-order 
tensor reads a’ and AT, respectively. The principal transposition 
of a fourth-order tensor is denoted via AT“ and the left and right 
transpositions are AT! and A'®. The linear mappings induced by second- 
order and fourth-order tensors are written as a = Cb and A = C[Bj, 
respectively. The composition of two second-order or two fourth-order 


tensors is denoted by AB and AB. The Frobenius inner product 
is denoted by A-B=tr(AB'). The tensor product operator & is 
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defined as (A® B) [C] =(B-C)A. Its symmetrized version ®s is 
defined via a &s b=(a®b+b8a)/2. We introduce the abbrevia- 
tion a®” = a ®a--- Qa (n repetitions). The material time derivative of 
a quantity w is expressed as ù = Dw/ Dt. 

We denote by Sym(d) the space of symmetric second-order tensors on Rë. 
The unit sphere in R? reads $?. The vector space of fourth-order tensors 
with minor symmetries (A = A™, A = A'®) is written as L(Sym(d)), 
whereas Sym(Sym(d)) denotes those fourth-order tensors that have 
minor and major symmetries (A = ATM), Generally, we use a {€1, €2, e3} 
Cartesian coordinate system on the (local) microscale and a {ez, €y, €z} 
Cartesian coordinate system on the (effective) macroscale. Details on 
further spaces of interest, domains of definition and corresponding 
explicit expressions are given upon their first appearance. 


Acronyms 

pCT Micro-computed tomography 
CDF Cumulative density function 
CDI Clausius-Duhem inequality 
CG Conjugate gradient 

Co Continuous 

CoDiCo Continuous-discontinuous 
DiCo Discontinuous 

FE Finite element 

FFT Fast Fourier-transform 

FRP Fiber reinforced polymer 
FRTS Fiber reinforced thermoset 
GRK Graduiertenkolleg 

GSM Generalized standard material 
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IRTG 
KKT 
LFT 
ODF 
OT 
PDF 
SMC 
UD 
UPPH 
YMS 


Latin letters 


a,b, A, B, 
a,b,c,... 
A,B,C,... 
A,B,C, 
c 

d 

f 

h 

m 

P 

q 

w 

z 

E 

G 

H 


International Research Training Group 
Karush-Kuhn-Tucker 

Long fiber reinforced thermoplastic 
Orientation distribution function 
Orientation tensor 

Probability density function 

Sheet molding compound 

Unidirectional 

Unsaturated polyester polyurethane hybrid 
Young’s modulus surface 


Scalar quantities 

Vectors 

Second-order tensors 

Fourth-order tensors 

Volume fraction 

Dimension 

Reduced damage-activation function 
Energy (density) accounting for degradation 
Power-law exponent 

Bayesian parameters 

Damage variable 

Free energy (density) 

Internal variables 

Young’s modulus 

Shear modulus 

Hardening parameter 
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Orthonormal basis vector 


m ao 


‚m,n,p Direction vectors 

Stress traction vector 
Displacement 

Current position 

Second-order orientation tensor 
Deformation gradient 
Displacement gradient 
Second-order identity tensor 


> sb hyp 8 e S 


Fourth-order orientation tensor 


Extraction tensor 


Q 


Stiffness tensor 


Žo Eshelby’s tensor 
Fourth-order identity tensor 


Po Polarization tensor 


P1, P2 Isotropic projectors 
S Compliance tensor 


Driving force for compliance evolution 
Set of active damage functions 
Expectation value 
Set of Gaussian normal distributions 
Set of damage variables 
Unit sphere 
(d) Set of symmetric tensors 
Set of sym. pos. definite fourth-order tensors 
Set of real numbers 


NP OLSRORS 


Set of internal variables 
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Greek letters 


an 


Driving force for damage evolution 
Scalar strain 

Consistency parameter 

Scalar stress 

Damage-activation threshold 
Backtracking factor 

Poisson’s ratio 

Orientation distribution function 
Damage-activation function 
Force potential 

Infinitesimal strain tensor 
Cauchy stress tensor 


Sub- and superscripts 


B’ Or ’ Om 

Lo ()r 

x? oF ’ (Ds 

I? On ’ (in ’ (div 


Bundle, Fiber, Matrix 
Longitudinal, Transverse 
Related to directions ez, €y, e; (orthotropy) 
Puck-type cases 

Parallel 

Perpendicular 

Related to normal stress 
Related to shear stress 
Symmetric 

Spherical 

Deviatoric 

Material time derivative 
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Operators 
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Composition of two second-order tensors 
Composition of tow fourth-order tensors 
Linear mapping by a second-order tensor 
Linear mapping by a fourth-order tensor 
Frobenius inner product of two tensors 
Tensor product of two vectors 

Symmetric tensor product 

Tensor product repeated n times 

Cross product of two vectors 

Transposition of a vector or second-order tensor 
Major transposition of a fourth-order tensor 
Left and right minor transpositions of a fourth- 
order tensor 

Inverse of a symmetric tensor 

Volume or ensemble average w. r. t. phase a 
Macaulay bracket 

Frobenius norm 

Eulerian divergence 

Eulerian gradient 

Symmetric part 

Skew-symmetric part 

Trace 

Diagonal tensor 


Chapter 2 


Fundamentals of continuum 
mechanics 


2.1 Kinematics 


In this chapter, we give an overview of continuum mechanical fun- 
damentals that provide the framework for our subsequent modeling 
approaches. We introduce kinematics and motivate the definition of 
displacement and infinitesimal strain. Furthermore, we provide an 
insight on the general formulation of balance equations and discuss 
specific, required balance equations. For further information and more 
detailed explanations, the reader is referred to the works of Truesdell 
and Toupin (1960), Holzapfel (2000), Haupt (2002) and Bertram (2011). 
A three-dimensional body, in the context of Cauchy continua, is defined 
as a set of material points that have three translational degrees of 
freedom each. The current position of a material point in such a body 
with volume v at time t is given by the position vector 


x= x(X,t), (2.1) 


where X denotes the position vector of the reference placement of the 
body at time t = tọ. The deformation x describes the change of the 
body due to boundary conditions, i.e., it describes the current shape 
of the body w.r.t. the corresponding reference placement. The inverse 
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mapping X = x~' (x,t) allows for the identification of the reference 
placement. The displacement (field) u (X, t) is the difference between 
the current and reference placement of all material points in a considered 
body 

u(X,t) =x (X,t)—- X. (2.2) 


A field quantity A can be defined in terms of the reference placement 
or the current placement, being called the Lagrangian A, (X,t) or the 
Eulerian Ag(a,t) description, respectively. The conversion between 


these two descriptions is performed via the deformation x and x”! 


Ar (x,t) = Ar lv” (æ,t),t), (2.3) 
Ay (X,t) = Ag (x(X,t),t). (2.4) 


The material time derivate of a field quantity (Truesdell and Toupin, 
1960) can be expressed via 
, OA, (X,t) _ OAg (a, t) OAg (ax, t) 


A= OE = ƏL + Ər Ë (2.5) 


with the velocity « = v(x,t) of the material point. Note that defini- 
tion (2.5) holds for any tensorial rank of the field quantity A in general. 
The deformation gradient F and the displacement gradient H, being 
spatial derivatives of the current position æ and the displacement u in a 
Lagrangian description, are defined as 


F (X,t) = Grad (z (X,t)) = Ox) (2.6) 
H (X,t) = Grad (u(X,t)) = mu) =F-I. (2.7) 


Despite its designation, the deformation gradient F also accounts for 
rigid body motions, even though the body does not deform. There exist 
several deformation measures, depending on the field of application 
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(Bertram, 2011). One deformation measure, taking into account the 
change of length and angles of line elements in a body, is Green’s strain 
tensor 
p=! (F'F = I) =. (H +H" + H"H) . (2.8) 
2 2 
The linearized version of the deformation measure is the (infinitesimal) 
strain tensor 
e (X,t) = sym(H) = ; (H J a ! (2.9) 


which is valid in the context of small deformations, i. e., the condition 


|| || = \/tr(HH") <1 (2.10) 


must hold (Bertram, 2011). The corresponding skew-symmetric part 
of the displacement gradient is the (infinitesimal) rotation tensor 
w = (H — H')/2=skw(H). 

Every second-order tensor can be unambiguously decomposed into a 
spherical and deviatoric part. For the strain tensor, we have 


e=e’+E, (2.11) 


where the spherical strain e° describes the volumetric change of the 
considered body (dilatation) and the deviatoric strain e’ gives the volume 
preserving change of shape of the body (distortion). We can extract 
these strains via two isotropic fourth-order projectors Pı = (I & I) /3 
and P» = IS — P4 as 


tr(e) tr(e) 


e° =P, [e] = I, e =Polel=e- I. (2.12) 


The fourth-order identity on symmetric second-order tensors is denoted 
by Is. 
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2.2 Balance equations 


2.2.1 General formulation 


A correct prediction of the (mechanical) behavior of a body is based 
on a correct application and evaluation of associated balance equations. 
To provide an overview of the most important equations, we introduce 
the general form of balance equations and subsequently derive specific 
applications or modifications, respectively. Balance equations associate 
the change of a field quantity A in a considered material volume v to the 
supply s and production p of this quantity, as well as the corresponding 
non-convective flux q across the boundary ðv (Truesdell and Toupin, 


1960), 
< [ Aew= fots) w+ f q: da. (2.13) 
dt v Ov 


v 
For the surface integral, the relation da = nda holds, with da being 
a surface element of the boundary ðv and n the outer surface normal 
unit vector. We can transform the integral form into a local form in 
regular points using the divergence theorem and Reynolds transport 
theorem (Bertram, 2011) 


^ + div (Av) = p + s + div (q). (2.14) 


2.2.2 Specific formulations 


Mass 


We obtain the balance of mass via the density A = ọ as 


0+ odiv (v) = 0, (2.15) 
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where the production p, supply s, and flux q are zero. We can calculate 
an alteration of mass via tr(e) = (dv — dvo) / dvo or o = oo (1 — tr(e)). 
In the framework of small deformations, i.e., condition (2.10) or 
rather ||e|| < 1 holds, the relation 9 = go is valid and all displacement 
fields automatically fulfill the balance of mass. 


Linear momentum 


The balance of linear momentum associates the rate of change of the 
density of linear momentum A = ov to the supply via volume forces ob 
and the flux in terms of surfaces forces t, leading to an integral form 


d 
Z| æw= f obav | tda. (2.16) 
dt v v ðv 


The production is zero. Using the Lemma of Cauchy t = øn, we can 
write the local form as 


ov = ob + div (ø). (2.17) 


For a quasi-static setting, the acceleration ù vanishes. 


Angular momentum 


Balancing the field quantity A = æ x ov yields the balance of angular 


momentum 


[ex ww [wx obav+ | x x tda, (2.18) 
v v Ov 


with the supply x x ob and the flux æ x t. For a Boltzmann continuum 
with vanishing moment densities, the quasi-static balance of angular 
moment in its local form yields the symmetry of the stress tensor 
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(Truesdell and Toupin, 1960) 
o=o". (2.19) 


Remark: The second-order Cauchy stress tensor ø describes the local 
stress state in a material point. Similar to the strain tensor, we can 
decompose the stress into a spherical part o° and a deviatoric part o’ 


oe=o’+c0, o° =Pilo], o =Po[o}. (2.20) 


Energy 


Relating the rate of change of total energy to the power furnished 
by volume ob- v and surface t- v forces, as well as volumetric heat 
supply os and heat flow across the surface q - n yields the integral form 
for the total energy 


zS jo A w= f (eb-v+ os) au+ | (t-v—q-n) da, 
dt v 2 v ðv 
(2.21) 


which is the first law of thermodynamics. The associated local form of 
the total energy in regular points is 


1 : 
oè + 50 (v-v) = ob- v + os + div (o'v) — div (q). (2.22) 


Note that the total energy is the superposition of internal energy and 
kinetic energy. We can derive the local form, i. e., the balance of internal 
energy, by taking into account the balance of mass and balance of linear 
momentum 

pe = os — div (q) +0 È. (2.23) 


The internal energy changes, due to the production of kinetic energy 
(density) ø - è due to stress power, heat flux across the surface div (q) 
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and heat supply in the volume os. In a similar way, we can obtain the 
balance of kinetic energy in its integral form via a multiplication of the 
balance of linear momentum with the velocity v (Iruesdell and Toupin, 
1960) and finally end up with the local form 


lolo o) = ob: v + div (o'v) —oa-€. (2.24) 


The flux of kinetic energy over the surface is given by the power of 
external forces div (a'v) = div (ov) and the supply of kinetic energy 
is ob - v. Energy in terms of stress power ø : € is exchanged between the 
internal energy and the kinetic energy. 


Entropy 


A rate of change of entropy n is evoked via entropy production opņ, 
entropy supply os, and entropy flux q,, -n 


d 
a fende= | ols +m) dw- f an nda, (2.25) 


Following the second law of thermodynamics, the entropy production is 
non-negative [, op, dv > 0 for any process. The local form of the entropy 
balance is given as (Coleman and Noll, 1963) 


on = 05n + OP — div (an) ‘ (2.26) 


Generally, the entropy supply and entropy flux are assumed to 
be s, = s/ and q, = q/0, see Coleman and Noll (1963). In terms of 
the local formulation, the second law of thermodynamics is given is 


py È 0. (2.27) 
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Utilizing the relation of the internal energy e, the entropy 7 and the 
Helmholtz free energy density w 


e=On+u, (2.28) 


as well as a proper Legendre transformation (Coleman and Gurtin, 1967), 
the local from the balance of internal energy and the local balance of 
entropy, we can derive the Clausius-Duhem inequality 


: 1 
o -È — ew — obn — 94° grad (0) > 0. (2.29) 
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Chapter 3 


Basics of SMC composites 


3.1 Manufacturing process 


Typically, SMC composites are manufactured via a compression molding 
process as shown in Fig. 3.1. Microstructure characteristics such as 
orientation or volume fraction are influenced by the manufacturing 
process leading to variations of the structural performance of produced 
parts. We give a brief overview of the manufacturing process established 
within IRTG 2078 to grasp the key aspects that affect the local and 
effective behavior of our considered SMC composites. 


Transfer of matrix to SMC line 


} CoFRTS bonded 
j+ with DICOFRTS 


CoDiCoFRTS 


+ DiCoFRTS with CoFRTS 


CoFRTS 


‘ea == == 


DiCoFRTS 


Compression Molding 


Figure 3.1: Manufacturing process of SMC composites (kindly provided by T. D. Pallicity 
(Karlsruhe Institute of Technology (KIT) - Institute of Engineering Mechanics)) 
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Depending on the intended application, the resin might be mixed with 
additives and fillers in a preliminary step, see the left hand side in Fig. 3.1. 
For the material at hand, we use an unsaturated polyester polyurethane 
hybrid (UPPH) resin with no fillers that possesses a stable bi-stage state 
allowing for a co-molding of glass and carbon fibers (Bücheler et al., 
2017). The so produced SMC composite is specifically designed to be 
used for the manufacturing of structural, load-bearing components. 
Subsequently, the mixed matrix is transfered to the SMC line and filled 
into two doctor boxes, see the middle in Fig. 3.1. These two doctor 
boxed spread the matrix on two different films, one upper and one lower 
film. Endless fiber rovings are cut to lengths of 25.4 mm (1 inch) in the 
chopper and are randomly dispersed on the lower film with matrix. 
Thereupon, the upper film with matrix is put on top of it, so that the 
fibers are encased in matrix. The produced DiCo pre-preg is fed into 
a carrier belt that induces relative motion between fibers and matrix 
and hence improves impregnation. Finally, the pre-preg is rolled up for 
better handling and storage. Partial curing (pre-curing) of the matrix can 
be induced by setting specific temperature conditions during storage. 
For the compression molding process of an SMC composite part, DiCo 
pre-pregs are cut and stacked to a specific size and thickness. This 
initial charge is placed on the tempered mold (typically 120° to 160°) 
at a certain position and typically does not fully cover the complete 
mold. Conceivably, Co patches are added at specific positions within a 
co-molding process, see the right hand side in Fig. 3.1. The closing of 
the press induces a plug-like flow of the DiCo SMC which ultimately 
fills the complete mold. The flow influences the fiber positions and 
orientations of the DiCo material and hence yields a process-dependent 
inhomogeneous orientation distribution in the manufactured part. If 
designed correctly, the Co patches remain at their assigned positions. 
Usually, it takes about 2 min. to 4 min. for the matrix to cure within 
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the mold. After that, the mold is opened and the manufactured part is 
removed via ejector pins. 


3.2 Fiber bundle structure’ 


Each of the utilized rovings in the manufacturing process (see Fig. 3.1) 
consists of about 225 aligned fibers. Despite being cut and dispersed on 
the film, these 225 fibers remain as bundles or tows, respectively. When 
encased in matrix and rolled up, these bundles remain intact, also during 
storage. Within the compression molding manufacturing process, the 
bundles are exposed to flow conditions in the mold and hence relocate 
and change their orientations. Nonetheless, due to various conditions, 
e. g., low viscosity of the matrix, plug-like flow behavior, outer (hot) 
lubrication layers and low shear forces, the fiber bundles mostly remain 
intact during the flow process. Therefore, typical SMC composite parts, 
manufactured via the discussed compression molding process, show a 
characteristic three-scale structure as presented in Fig. 1.6. 


(a) Volumetric „CT image (b) Identified fiber bundles (c) Generated microstructure 
(Schöttl et al., 2021b) (Schöttl et al., 2021b) (Görthofer et al., 2020) 


Figure 3.2: Real fiber bundle microstructure of an SMC composite (Schöttl et al., 2020; 
Schöttl et al., 2021b) and generated microstructure (Görthofer et al., 2020) 


1 This section is based on excerpts of the publication “A computational multiscale model 
for anisotropic failure of SMC composites” (Görthofer et al., 2022a) 
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Dedicated algorithms (Pinter et al., 2018; Schöttl et al., 2020; Schottl et al., 
2021b) reveal the fiber bundle structure of a manufactured component 
on pCT images using tracking and clustering techniques, see Fig. 3.2. 
Similar microstructure characteristics are observed for mesostructure 
analysis via optical microscopy on polished SMC composite sample 
surfaces (Chen et al., 2018b; Li et al., 2018a) and accounted for in 
hierarchical approaches (Anagnostou et al., 2018; Fitoussi et al., 2005). 


3.3 Orientation tensors’ 


To take the orientation state of our considered SMC composite parts into 
account within the material modeling approaches, we need to depict the 
orientation via mathematical descriptions. For a volume containing 
identically shaped cylindrical fibers, following Kanatani (1984) and 
Advani and Tucker (1987), we describe the fiber orientation state by 
an orientation distribution function (ODF) y : S”71 — R, defined on the 
unit sphere S”~' in R” (n = 2,3). The fraction of fibers pointing in a set 
of directions D C S”~' is determined by 


sn HP) ds, 6-1) 


where ds denotes the surface element of the (n — 1)-sphere, and |S”~'| 
denotes the surface measure of S”~'. Considering n = 3 dimensions, 
we have a surface element of the unit sphere ds = sin 0 dé dy in terms of 
spherical coordinates {r, 0, y}. For physical reasons, the ODF is assumed 
non-negative, normalized (i.e. integrates to unity) and symmetric, i.e., 
w(—p) = v(p) holds for any p € S”~!. 


1 This section is based on excerpts of the publication “Computational homogenization of 
SMC composites based on high fidelity representative unit cells” (Görthofer et al., 2020) 
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The particular forms of Eq. (3.1) for the two- and three-dimensional case 
(n = 2,3) are given by 


1 27 


»b(p) dp, p= (cosp,siny,0)', and (3.2) 
y=0 


1 27 T 
/ U(p)sin@ dédy, p = (sin 9 cos ọ, sin @ sin Y, cos a)". 
yg 


(3.3) 


Hereby, the planar and spatial unit spheres and their surface measures 
are defined via the radius r. The corresponding ODFs ı(p) need to be 
planar or spatial distribution functions, respectively. 

Based on the ODE, orientation tensors (OTs) may be defined by forming 
dyadic products of the directions p and a subsequent integration over all 
directions in the (n — 1)-sphere. Due to symmetry of the ODF, the OTs 
of odd-order vanish. The even-order integrals yield totally symmetric 
and positive definite tensors with trace 1. The most commonly used OTs 
are the second- and fourth-order OTs A € R"*" and A € Rn 


A= a OPEP, (3.4) 


A= Yp) pS ppp ds, (3.5) 
Sn-1 
where we suppress the dependence on n for notational clarity. A general 
second-order OT in three spatial dimensions (defined by an orthonormal 
basis {eı, €2, e3} may be represented as 


Aı Aj. Ais 
A= |Aıa Ago Az (3.6) 
Aı3 Ao3 Az33 
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where the constraint 
tr(A) = Aıı + A22 + A33 = 1 (3.7) 


is a consequence of fgn-ı Y(p) ds = 1. 

Due to the compression molding manufacturing process, the fiber 
orientation distribution turns out to be almost perfectly planar for SMC 
composites, i.e., the components pointing in the normal direction of 
the produced part vanish. Thus, we may describe such an orientation 
distribution by second-order OTs on the unit circle S?. Furthermore, we 
can diagonalize the second-order OT A by a proper rotation matrix R 


Ai Aıa 
Aıa Ago 


A 0 
0 A 


R' (3.8) 


with eigenvalues A; and As of the planar second-order OT A. 


3.4 The exact closure approximation in two spa- 
tial dimensions’ 


The principal advantage of fiber orientation tensors compared to 
working with the full ODF is the reduced complexity (and, thus, the 
reduced number of degrees of freedom) necessary for numerical process 
simulations. Unfortunately, the second-order fiber orientation tensor 
does not determine the effective elastic properties of a fiber reinforced 
composite, in general. On the other hand, the phase space of fourth-order 
fiber orientation tensors is extremely complex, cf. Linn (2005). As a 
work-around, Tucker and co-workers (Advani and Tucker, 1987; Cintra 
and Tucker, 1995) pioneered working with closure approximations, i. e., 


1 This section is based on excerpts of the publication “Computational homogenization of 
SMC composites based on high fidelity representative unit cells” (Görthofer et al., 2020) 
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functions g mapping second-order OTs to (reasonable) fourth-order OTs, 
A = g( A). The first closure approximations, like the linear, quadratic 
and hybrid closure approximations (Advani and Tucker, 1987) violated 
physical principles, which comes from the fact that these closure 
approximations could generate fourth-order OTs that do not derive 
from an ODF, i. e., there was no 7 : S"~! + R, s. t. Eq. (3.5) holds. This 
problem could be resolved by factoring closure approximations, i. e., to 
seek mappings h associating an ODF y to a second-order OT A. The 
resulting closure approximation subsequently computes as 


A= = U(p)p®pBp@pads, where w(p) =h(A). (3.9) 


Thus, finding closure approximations for OTs reduces to seeking a 
suitable class of ODFs parameterized by second-order OTs. 

As shown by Verleye and Dupret (1993), the family of angular central 
Gaussian distribution functions (Tyler, 1987) 


1 ñ 
YB(p) = pmj -p@p)?, pes", (3.10) 


parametrized by a positive definite symmetric matrix B € R"*” with 
determinant 1, constitute exact solutions for the ensemble Jeffery’s 
equation (Jeffery, 1923), and thus, to the Folgar-Tucker equation (Folgar 
and Tucker, 1984) with vanishing Folgar-Tucker diffusivity, which 
model the fiber orientation dynamics during injection and compression 
molding. 

For n = 3, Montgomery-Smith et al. (2011), extending earlier work 
of Verleye and Dupret (1993), noticed that there is a one to one 
correspondence of admissible B’s to second-order OTs A, via the 
association derived from Eq. (3.4) with a distribution function Y = we 


A= N Yg(p) p® p ds. (3.11) 
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Unfortunately, determining B for given A involves solving a system of 
elliptic integrals which becomes ill-conditioned for det(A) < 1. 
In this work, we are concerned with the two-dimensional case, i.e., 


n = 2, and show that there are simple explicit formulae for B in terms 
of A (and vice versa). Also, the fourth-order OT as defined in Eq. (3.5) 


a= | vB(p)p® p® p®p ds (3.12) 
sr—-1 


can be computed explicitly. To be more precise, the following holds: Let 
B be a diagonal matrix B = diag( B1, B2) with Bı > 0 and By = 1/B.. 
Then, A computed from (3.11) is a diagonal matrix A = diag( A1, A2) 


with entries 1 


A en, 
1 1+ Bə 


and A= 


1 
STB (3.13) 


Furthermore, the (independent) components of the fourth-order OT 
(A = Aijkl €i D €j Q ekx 8 e,) as defined in Eq. (3.12) read 


1 
Ann = z410 + Aj), A2222 = A, ) 
(8.14 


il 
A1122 = z442, A1112 = 0, Aj222 = 0. 


Hereby, the components are invariant under permutations of the indices. 
A derivation is contained in Appendix A. 

Using a proper rotation, these formulae for the exact closure can also be 
used for non-diagonal OTs. Another consequence is the identity 


(3.15) 
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which only holds for n < 2. To establish this identity, consider the 
diagonal case. Then, the latter identity is equivalent to 


1 1 
a =A, and ao = Ap. (3.16) 
VB, VB2 VBı VB 
However, using Bı B2 = 1, we deduce 
1 
VBı vV Bə vV Bı Bo 1 


(3.17) 


gat Te BB Bti8B Bisel’ 


which equals Aı by the first equation in (3.13). The second identity is 
shown analogously. 

A most widely used closure approximation is the quadratic approx- 
imation (Doi, 1981; Advani and Tucker, 1990). Hereby, the fourth- 
order orientation tensor A is approximated via the dyadic product 
of the second-order orientation tensor A with itself (A = A® A). 
With A = diag(Aı, A2), the entries of A are 


2 2 
A1111 = Aj, A2222 = A3, 


(3.18) 
A1122 = A; Ao, A1112 = 0, A1222 = 0. 


When we compare the two closure approximations, we see the superi- 
ority of the exact closure Eq. (3.14) over the quadratic closure Eq. (3.18). 
The component A1122, €. g., is off by a factor of 2 and therefore not 
predicted appropriately via the quadratic closure. 
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Chapter 4 


Computational homogenization 
of SMC composites based on high 
fidelity representative unit cells’ 


4.1 Introduction 


4.1.1 Research contributions 


Sheet molding compound (SMC) composites combine high lightweight 
potential with excellent formability and are frequently used in industrial 
applications. To reduce safety factors in dimensioning SMC parts, 
the influence of processing parameters and stochastic variation of mi- 
crostructural and physical properties needs to be quantified accurately. 
Taking into account the inherent three-scale structure of SMC, we 
improve the microstructure generator of Chen et al. (2018b) in various 
respects. Firstly, we consistently rely upon state-of-the-art closure 
approximations for the fourth-order fiber orientation tensor. More 
precisely, we show that for planar fiber orientation state, there is an 
explicit formula for the fast exact closure approximation of Montgomery- 
Smith et al. (2011). Secondly, we exploit the use of quasi-random 
numbers in sampling the fiber orientation distribution, leading to 


1 This chapter is based on the publication “Computational homogenization of SMC 
composites based on high fidelity representative unit cells” (Görthofer et al., 2020) 
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dramatic improvements in accuracy compared to pseudo-random Monte 
Carlo sampling. Last but not least, we rely upon fast Fourier transform 
based methods for rapid computational homogenization. 

With these methodological improvements at hand, we thoroughly inves- 
tigate the influence of the mechanical and microstructural parameters 
on the effective elastic properties of SMC composites, and compare the 
results to direct numerical simulations on large scale digital volume 
images and mean-field estimates. 


4.1.2 Chapter structure 


In this chapter, we are concerned with E-glass fiber reinforced SMC 
based on an unsaturated polyester polyurethane hybrid (UPPH) resin 
without fillers. This material class has been investigated in terms of its 
molding capabilities (Hohberg et al., 2017) and the resulting mechanical 
properties, in particular the viscoelastic characteristics (Kehrer et al., 
2018), puncture properties (Trauth et al., 2019), damaging behavior 
(Trauth et al., 2018; Schemmann et al., 2018b), tension-compression 
asymmetry and bi-axial characterization (Schemmann et al., 2018c). 
Due to the high bundle aspect ratio and the strong manufacturing 
process dependence of the fiber orientation, the stochastic bundle 
orientation selection strategy of Chen et al. (2018b) proved insufficient 
for our purposes. More precisely, we improve the bundle deposition 
algorithm in various regards. 

First, we consistently rely upon the powerful exact closure approxi- 
mation of Montgomery-Smith et al. (2011). We derive a closed-form 
expression for the exact closure of planar fiber (bundle) orientations. This 
result, described in Sec. 3.4 and Appendix A, is surprising, because for 
general three dimensional fiber orientation states, a system of equations 
involving non-standard elliptic integrals needs to be solved numerically. 
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Our second improvement concerns the strategy for choosing the fiber 
bundle directions. In contrast to the typical pseudo-random selection 
strategy, for instance relied upon by Chen et al. (2018b), we advocate 
using quasi-random sampling in the form of Sobol’s sequence (Sobol’, 
1967). The improvements in accuracy are dramatic, cf. Sec. 4.2.3, enabling 
us to use only a few bundles to gain high fidelity in targeted fiber 
(bundle) orientation tensors. 

Last, but not least, we consistently rely upon modern FFT-based 
computational techniques (Moulinec and Suquet, 1994; 1998; Schneider, 
2021), permitting us to quickly and thoroughly investigate the sensitivity 
of the effective elastic properties w. r. t. the multitude of both mechanical 
and geometrical parameters entering the SMC unit cell generator 
based computational multiscale strategy, cf. Sec. 4.4. Of course, we 
pay close attention to working with representative results, both in 
terms of spatial resolution and representativity of the underlying 
volume elements. Furthermore, we compare our findings to mean-field 
predictions and full-field simulation results on CT image data of an 
SMC microstructure. 


4.2 Generating SMC unit cells 


4.2.1 SMC materials and processing 


SMC composites are manufactured via compression molding process, cf. 
Bücheler et al. (2017) and Bücheler (2018). Initially, fiber tows consisting 
of approximately 225 individual, aligned continuous fibers with a diam- 
eter of approximately 13.5 um are cut to a length of 1 inch = 25.4 mm. 
These fiber bundles are randomly dispensed onto a matrix resin foil, 
enveloped by another matrix resin foil, rolled up onto sheets and 
pre-cured. From these pre-cured sheets initial charges can be cut to 
be used within the compression molding process. The measured fiber 
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volume fraction cr we consider for this work is 25 %, cf. Görthofer et al. 
(2019b). 


(a) Cross-section of SMC composite unit (b) Bundle placement in a voxelized unit cell 
cell showing the layerwise structure in z- 
direction 


Figure 4.1: Layerwise structure of SMC composite and algorithmic bundle placement to 
capture the layerwise structure 


During the compression molding manufacturing process, the matrix 
resin is heated, flows and fills the mold, and is compressed (Hohberg 
et al., 2016; 2017). Hereby, the individual tows of 225 fibers are carried 
by the flowing resin and are flattened but remain otherwise intact. The 
resulting SMC microstructure, cf. Fig. 4.1a, is composed of closely packed 
fiber “bundles”. These bundles remain almost unbent, hence we model 
them as planar. The fibers within a bundle are almost perfectly aligned. 
Thus, instead of describing the fiber orientation state by fiber OTs, we 
may also rely upon a description in terms of “bundle OTs” instead. 


4.2.2 The successive bundle deposition method 


As the point of departure for generating representative SMC unit cells 
serves the successive bundle deposition method described in Chen et al. 
(2018b), originally developed for carbon fiber SMC. Our method differs 
in the way the orientations are sampled, cf. Sec. 4.2.3. 
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In the algorithm, cuboids serve as the basic geometric primitives for 
modeling the fiber bundles, described by a bundle length lg, bundle 
width wg and bundle height hg (which we assume identical for the 
complete microstructure). These bundles are deposited by a random 
sequential addition (RSA) technique (Feder, 1980) into a square unit cell 
with dimension L, = Ly = L and periodic boundary conditions. For 
physical reasons, the bundles are prohibited to overlap. To reach the 
fiber volume content typically encountered in industrial applications (up 
to 60 % bundle volume fraction), bundles are allowed to be displaced 
into the next layer (but not more than a single layer). This is shown 
in Fig. 4.1b, which represents an actual result of the bundle deposition 
algorithm. For efficient collision handling, a background voxel grid with 
voxel spacing v is superimposed. 


input: 


unit cell: Ly, Ly Niayer V 


bundles fh, cat fhe * random bundle midpoint 


* quasi-random direction q 


layers * calculate bundle direction 
p based on q and OTA 
for each layer from 1 to Nıayer: | 


+ read layer information 


* while ny < MAX zies! 
m> * place bundle * voxelize bundle 

* check if feasible | + check if bundle 
can be placed in 
current layer 
and/or layer above 


e while current c < cp: 


* generate new bundle based 
on OT A of current layer 


yes 


aS 


done 


« place bundle (until it fits) an 


e voxelize bundle and calculate 
new current c and A 


+ change midpoint 


v 


| * voxelize bundle 3 
+ place bundle in curren We 
output: layer and/or layer above 
voxelized unit cell V * save orientation 


* microstructure 
* orientation 


Figure 4.2: Schematic flowchart of the key steps of the algorithm 
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Furthermore, the number of layers and a bundle volume fraction cg plus 
a bundle OT of second-order A can be defined for each layer. The result 
of the algorithm is a binary description of the microstructure on a voxel 
grid, together with bundle directions for the voxels containing bundles. 
The key steps of the algorithm are the generation of the bundles 
as abstract objects, determination of the directions according to the 
prescribed OTs, placement of the bundles in the corresponding layers 
and the voxelization of the bundles in the unit cell. The algorithm 
proceeds layer by layer, generating bundles sequentially, each of which 
samples random midpoints until the bundle does not overlap with 
the already existing bundles. Here, overlapping is defined taking into 
account a possible displacement into the next layer. Figure 4.2 presents 
a flow chart of the algorithm. 


(a) Ay = 1.0 (b) Ay = 0.5 (c) Ay = 0.2 (d) Aı = 0.0 


Figure 4.3: Generated unit cells with different orientations 


Examples of generated unit cells based on the presented algorithm are 
shown in Fig. 4.3. For better illustration of the placed bundles, the 
number of shown layers is limited. Blue bundles are placed into a 
white matrix, and the shadings of blue indicate the different layers. A 
bundle volume fraction of cg = 50% was chosen for all layers. The 
desired volume fraction is matched with a relative error of e < 0.5%. As 
error measurement we normalize the difference of the desired and the 
actual volume fraction w. r. t. the desired volume fraction. The presented 
microstructures vary in the desired orientation distribution going from 
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unidirectional in x-direction (Fig. 4.3a) via planar isotropic in the z-plane 
(Fig. 4.3b) to unidirectional in y-direction (Fig. 4.3d) for all layers. Using 
a similar error measurement as for the volume fraction, the orientations 
are matched with a relative error of e < 1.5%. 


4.2.3 Quasi-random orientation sampling 


The key addition to the microstructure generation method of Chen et al. 
(2018b) is the way in which orientations are sampled. There are two 
requirements. On the one hand, for a given planar OT A, the drawn 
directions p4, ..., pg should match the fourth-order OT A determined 
by the exact closure to high precision, i. e., the difference 


K 
CESSAT ATATA > min. (4.1) 
i=1 
should be as small as possible. Hereby, we measure the difference via the 
Frobenian norm. Secondly, the first property should hold even though 
the number of drawn directions K is not known a priori. The latter 
results from the bundle deposition method described in Sec. 4.2.2 - as 
parts of each bundle may be displaced to the next layer, the number of 
bundles to be drawn to match the desired fiber bundle content has to be 
determined dynamically during the algorithm’s performance. 
Notice that if we knew K beforehand we could draw the direc- 
tions p,,...,P, randomly and perform an optimization scheme to make 
(4.1) as small as desired. 
Goldberg et al. (2017) noticed that the ODF determined by the exact clo- 
sure approximation arises as a transformation of a uniformly distributed 
random variable on the unit sphere S"~'. More precisely, suppose a 
second-order OT A with det(A) > 0 is given. Suppose that B is the uni- 
modular matrix realizing the angular central Gaussian distribution wg 
with second-order OT A. Then, if qis a multivariate random variable 
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of vectors uniformly distributed on the unit sphere S”"~', the random 
variable p = B~?q / IB":q| is wp distributed. Here, B~? is the inverse 
of the square root of the matrix B which may be computed, for instance, 
by a spectral decomposition. 

For the situation at hand, n = 2 holds, and the previous formula 
simplifies. Suppose p is a random variable uniformly distributed on the 
unit circle S$. Then, by the scaling identity, the random variable (3.15), 


Aq 


P= Tagl a 


is Wp-distributed. Thus, surprisingly, for the planar exact closure 
approximation we do not need to determine the matrix B to transform 
a uniform distribution to the ACG distribution. 
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Figure 4.4: Accuracy of generated orientation 
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With this insight at hand we would be able to generate bundle 
directions by sampling from a uniform distribution on the unit circle, 
or, equivalently, on the interval [0, 1]. Such a method would be similar 
to the method Chen et al. (2018b) use, with the additional benefit of 
using the exact closure of Montgomery-Smith et al. (2011). To improve 
the convergence rate of this Monte-Carlo-type sampling, we use the 
low-discrepancy sequence introduced by Sobol’ (1967) for drawing 
points in [0,1]. 

To illustrate the increase in accuracy, relative deviations based on the 
Frobenian norm between the desired second-order OT and the realized 
OT of a generated unit cell are shown in Fig. 4.4. In the diagrams, several 
desired orientations, represented by the largest eigenvalue A, of the 
second-order OT A, are investigated, together with a single stochastic 
realization. 

Quasi-random sampling leads to a more-or-less monotonic decrease 
of the error for increased number of bundles. In contrast, for pseudo- 
random sampling, the error may be strongly non-monotonic. Further- 
more, the reached accuracy for quasi-random sampling is much higher 
than for pseudo-random sampling, in accordance to theory. 


4.3 Computational setup 


4.3.1 Implementation, soft- and hardware 


The SMC unit cell generator (Sec. 4.2) was implemented in Python 2.7 
with Cython extensions. For the full-field simulations, we rely upon a 
fast Fourier-transform based computational micromechanics code, as 
described by Schneider (2018), using the staggered grid discretization 
(Schneider et al., 2016b) and the conjugate gradient method of Zeman 
et al. (2010). The computations ran on 6 threads on a desktop computer 
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with 32 GB RAM and an Intel i7-8700K CPU with 6 cores and a clock 
rate of 3.7 GHz. 


The maximum runtime for a voxel mesh size of 1500 x 1500 x 42 voxels 
(= 95 - 10° voxels) was about 6 - 5 min = 30 min for computing all com- 
ponents of the effective stiffness. The average run time on a voxel mesh 
with 750 x 750 x 21 voxels (~ 12 - 10° voxels) was about 6 - 30 s = 3 min. 
The mean runtime in order to generate a unit cell with a bundle volume 
fraction of cg = 50% was about 3s to 10s. 

Our approach is exemplified in Fig. 4.5, where both a generated SMC 
microstructure and the local von Mises equivalent strain and stress fields 
are shown for a uniaxial in-plane extensional experiment. Apparently, 
the fiber bundles pointing in loading direction carry the load, whereas 
the stress is much lower in the other bundles. Strain concentrations arise 
at the fiber bundle tips, but appear not to be very pronounced. 


40 MPa 
20 MPa 
0 MPa 


(a) Generated microstructure (b) Local von Mises strain field (c) Local von Mises stress field 


0.30 % 


0.15% 


0.00 % 


Figure 4.5: Generated SMC microstructure and computed local strain/stress fields for 0.1 % 
applied tensile strain in horizontal direction 
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4.3.2 Elastic properties of the constituent phases and the 
bundles 


The isotropic phase properties of the considered UPPH matrix resin 
system (Bticheler, 2018) and E-glass fibers are listed in Tab. 4.1. We 
computed the transversely isotropic fiber bundle properties via two 
different complementary methods. On the one hand we used numerical 
full-field homogenization of a representative bundle containing 225 
unidirectionally aligned fibers with a fiber volume content within 
a bundle of cr = 50%, see Sec. 4.2.1. A bundle volume fraction 
of cg = 50% in combination with cpg = 50% yields the corresponding 
fiber volume fraction of cr = 25% as presented in Sec. 4.2.1. On the 
other hand, we used an orientation averaged mean-field estimate of 
Mori-Tanaka type (Mori and Tanaka, 1973; Duschlbauer et al., 2003), 
as described in Sec. 2.1.4 of Brylka (2017). Both methods lead to 
approximately the same transversely isotropic properties of the fiber 
bundles. These computed elastic engineering constants, with "L" 
standing for longitudinal and "T" for transverse, are listed in Tab. 4.1. 


Table 4.1: Elastic material parameters for the considered SMC composite 


E-glass fibers UPPH matrix Fiber bundles 
(Kehrer et al., 2018) (Kehrer et al., 2018) 
Eriso = 72 GPa EM iso —=3.4GPa EBL = 37.73 GPa 
ET = 10.33 GPa 
VE iso = 0.22 UM, iso = 0.385 VB,TT = 0.477 
vg Lr = 0.292 


GEiso = 29.51 GPa GM iso = 1.23 GPa GET = 3.58 GPa 
GBLr = 3.64 GPa 
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4.3.3 On isotropic, transversely isotropic and orthotropic 
approximations of the effective elastic properties of 
the SMC composite 


Using fast Fourier transform based computational micromechanics 
enables computing the effective stiffness tensor for a generated SMC 
unit cell. Interpreting the 21 resulting constants may be cumbersome, 
in general. As we always align the axes of the unit cell with the 
principal axes of the OT, it may be tempting to rely upon an orthotropic 
approximation of the computed effective elastic tensor, and to interpret 
the associated engineering constants instead (see, e. g., Cowin (1985) or 
Bohlke (2001)). In this paragraph, we investigate numerically whether 
this is plausible, also including the isotropic and transversely isotropic 
approximation errors. 

As presented in Sec. 4.3.2, we use a fiber volume fraction of cr = 25% 
and a corresponding bundle volume fraction of cg = 50 %, as well as the 
phase and bundle properties listed in Tab. 4.1 for the considered SMC 
composite. The underlying unit cell parameters are listed in Tab. 4.2. A 
detailed and structured derivation of these parameters is presented in 
Sec. 4.3.4. 

For this study, we varied the A; component of the planar second-order 
OT A, as presented in Eq. (3.8). We implemented approximations for 
several classes of anisotropy, cf. Browaeys and Chevrot (2004). The 
resulting relative approximation errors of the stiffness matrices in Voigt- 
Mandel notation using the Frobenian norm depending on the orientation 
of the unit cell is shown in Fig. 4.6. Furthemore, representative plots 
of the Young’s modulus of the effective stiffness following Bohl ke and 
Brüggemann (2001) are included. 

Not surprisingly, an isotropic approximation of the elastic properties 
of an SMC composite is very inaccurate, leading to a minimum error 
of 17%. Instead, an orthotropic approximation of the SMC composite 
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Figure 4.6: Comparison of different approximations for the effective stiffness of SMC 
composites, including some exemplary Young’s modulus plots 


stiffness tensor leads to an approximation error of less than 1%, 
regardless of the orientation. For particular orientations, namely 
unidirectional orientation of all bundles in x-direction, in y-direction 
and a planar isotropic orientation in the z-plane, the effective stiffness 
of the unit cell is transversely isotropic. Depending on the orientation, 
the transversely isotropic axis switches from the y-axis to the z-axis to 
the x-axis. 

Throughout the following sections, we will analyze the transversely 
isotropic parameters EL, Er, Grr, Gir, vrr and vrr (Boehler, 1987), if 
possible. Hereby, "L" indicates the longitudinal and "T" the transverse 
direction. Otherwise, we will consider the orthotropic elastic parameters. 
We adapt the commonly used engineering notation Ex, Ey, Ey, Gyz, 
Gyz, Gxys Vyzs Vxzr Vays Vays Vex and Vyx (Boehler, 1987). Accordingly, we 
have 1 = x, 2 = y and 3 = zw.r.t. the introduced notation. For the sake 
of completeness, we use 6 transversely isotropic parameters and 12 
orthotropic parameters. 
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4.3.4 Determining the proper resolution, bundle aspect 
ratio as well as unit cell size and thickness 


Voxel mesh 


Unit cell 


Figure 4.7: Geometric parameters needed to ensure representativity of the unit cell, (a) 
unit cell dimensions, (b) bundle dimensions, (c) voxel mesh size 


To ensure reproducibility of the computational results, several sources 
of errors need to be taken care of, cf. Fig. 4.7. 


1. For fixed geometry, the resolution needs to be sufficiently fine to 
approximate the “continuous” results to sufficient precision. The 
geometry produced by the bundle deposition method of Sec. 4.2 is 
naturally given on a voxel grid. To study the influence of the resolu- 
tion, we introduce a positive integer M which we call amplification 
factor. For given M > 1, each original voxel is subdivided into M 3 
congruent voxels for the numerical investigations. 

2. For fixed microstructural parameters and sufficient resolution, the 
unit cell needs to be large enough to capture the statistical variation 
and the length correlation of the inclusions. Such unit cells are called 
representative volume elements (RVEs), cf. Gusev (1997). For the 
problem at hand, we need to investigate the dependence of effective 
elastic properties on the in-plane length L, = Ly = L and the number 
of layers Njayer (encoding the out-of-plane size). 


4.3 Computational setup 


3. For short fiber reinforced composites, it is known that the effective 
elastic properties saturate (Sun and Vaidya, 1996) for increasing 
aspect ratio, i.e., the quotient of fiber length and diameter. Thus, 
for computational purposes it is often sufficient to work with fibers 
that are much shorter than in reality. A natural strategy to determine 
this length is by successively increasing the fibers’ aspect ratio until 
the effective elastic properties do not change significantly any more. 


For the SMC microstructures at hand, we follow a similar strategy. 
From CT images of SMC (Le et al., 2008), such as Fig. 4.1a (see also 
Fig. 4.16a), we get a good estimate for the dimensions of a bundle’s 
cross section, i.e., the bundle width wg and height hg. However, 
the aspect ratio lg /wg, where lg denotes the bundle length, appears 
rather large. In general, the bundle aspect ratio of length to width 
varies approximately between lg/wg ~ 15 — 45, normalized w. r. t. the 
height hg. Thus, we investigate whether a smaller effective bundle 
aspect ratio is sufficient. 


Notice that effective elastic properties computed on a unit cell are 
independent of the unit cell length scale, i. e., there is no internal length 
scale present in the effective properties. Thus, the absolute length scale 
of the microstructure does not matter. Instead, relative length scales 
are of importance, i.e., the unit cell dimensions relative to the bundle 
dimensions, or the voxel size compared to the bundle size. For the 
investigation at hand, we choose the voxel size v initially produced by 
the bundle deposition method, cf. Fig. 4.2, as our basic length unit. 
Determining the various parameters satisfying all three stated require- 
ments involves an iterative process based, to some extent, on trial and 
error. For the purpose of presentation, we chose a set of “typical” 
parameters resulting from such an iterative process and show, by means 
of sensitivity studies, that their values ensure representative mechanical 
results (within engineering accuracy). 
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We consider a bundle volume fraction of cg = 50% as well as a planar 
isotropic orientation with A; = 0.5, representing an SMC microstructure 
frequently encountered in our application. To study the influence of 
the in-plane length L, the effective transversely isotropic engineering 
constants of an SMC composite are shown in Fig. 4.8, measured in 
the edge length of the original voxels, as discussed above. For the 
investigation, all remaining parameters were kept constant at Mayer = 7, 
lg = 50 voxels, wg = 5 voxels and hg = 1 voxel with a mesh resolution 
amplification factor of M = 3, using the phase properties given in 
Tab. 4.1. Fig. 4.8 contains both absolute values and relative errors, 
computed relative to the result computed at highest in-plane length 


considered. 
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Figure 4.8: Effective transversely isotropic engineering constants computed on unit cells 
with increasing in-plane length L for a planar isotropic orientation distribution 


Even for the lowest in-plane-length of 100 voxels, the engineering 
constants deviate less than 2 % from their “converged” counterparts. We 
see that an in-plane length of L = 150 voxels is sufficient for ensuring a 
relative error less than 0.5 % for all engineering constants and the chosen 
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geometric parameter set. Put differently, it suffices to work with a unit 
cell length which is three times the bundle length to obtain representative 
results. To ensure a degree of flexibility, we chose the unit cell in-plane 
length as L = 250 voxels. 

Similar studies were conducted to determine the necessary resolution 
and the magnitude of the remaining geometric parameters, the layer 
count Niayer and the bundle aspect ratio lg /ws, which ensure represen- 
tativity of the unit cells and the ensuing computational results. Fig. 4.9 
contains the corresponding graphs for increasing amplification factor, 
layer count and bundle aspect ratio. The determined parameter values 
ensuring representative effective elastic properties are listed in Tab. 4.2. 
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Figure 4.9: Relative error diagrams for increasing resolution, layer count and bundle aspect 
ratio 


In accordance with Tab. 4.2, for bundle dimensions to 50 x 5 x 1 voxels, 
a planar orientation state A; = 0.5 and a bundle volume fraction 
of cg = 50%, a relative error less than 1 % is ensured by choosing M > 3, 
cf. Fig. 4.9a. 

To minimize the influence of the unit cell dimensions on the computed 
effective stiffness, for studying the necessary bundle aspect ratio we 
choose the dimensions as L x L x Mayer = 500 x 500 x 7 voxels with 
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an amplification factor of M =3. The orientation remains planar 
isotropic Aı = 0.5 and the bundle volume fraction is cg = 50%. For 
an aspect ratio of 10 and higher, the relative error is less than 1% 
and the effective elastic parameters do not change significantly, as 
we can observe in Fig. 4.9c. Therefore we choose said aspect ratio 
of lg/wg = 10. Thus, bundles with smaller aspect ratio than for real SMC 
microstructures suffice. As discussed in Sec. 4.3.4, in SMC composites, 
the typical bundle aspect ratio range from approximately 15 to 45. 


Depending on the manufacturing process, the amount of processed 
material and the geometrical specifications of the component, the 
number of layers Mayer of a SMC composite varies. The pCT scan 
in Fig. 4.1a shows a number of layers Njayer of about 15 to 25. The 
question arises whether we need to take into account all layers in 
order to correctly predict the effective stiffness. For this last study, 
we use a unit cell in-plane length of L = 250 voxels, bundle dimensions 
of lg x wg x hg = 50 x 5 x 1 voxels and an amplification factor of M = 3. 
The orientation distribution is planar isotropic A; = 0.5 in each layer 
and the bundle volume fraction is set to cg = 50%. For Mayer = 7 layers 
or more, the effective elastic parameters barely change, cf. Fig. 4.9b, i. e., 
we ensure a relative error not exceeding 1%. 


Table 4.2: Geometric parameters chosen to ensure representativity of the generated unit 
cells and the typical parameters of the considered unit cells 


Unit cell dimensions Lx LX Mayer 250 x 250 x 7 voxels 
Bundle dimensions lg x wp x hp 50 x 5 x 1 voxels 
Amplification factor M 3 

Fiber volume fraction CE 25% 

Bundle volume fraction cg 50% 

Fiber fraction in bundle cpg 50% 

Orientation A=diag(A1,A2) Aı = A2 = 0.5 
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4.4 Multiscale simulation-based studies of ef- 
fective elastic properties 


4.4.1 Influence of volume fraction 


In this section we rely upon the material and geometric parameters fixed 
in Tab. 4.1 and Tab. 4.2. For the sensitivity analysis, we modify material 
parameters such as the orientation (by varying the A, component of 
the second-order OT A), the volume fraction cg, as well as the elastic 
properties of matrix Em, vm and fibers Er, vr. Comparing the resulting 
effective elastic parameters helps distinguishing important properties 
from unimportant factors. Due to the multitude of parameters entering, 
the latter is crucial for improving corresponding experimental setups, 
for instance. 
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Figure 4.10: Computed effective engineering constants depending on the volume fraction 
for cg € [40 %, 60 %] 


First, we computed the effective elastic parameters with varying volume 
fraction cg. We chose a planar isotropic orientation, i.e., A; = 0.5. 
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Thus, the effective stiffness tensor is transversely isotropic w.r.t. the 
z-direction. For volume fractions ranging from cg = 40% to cg = 60%, 
the corresponding effective engineering constants are plotted in Fig. 4.10. 
These computed stiffnesses parameters are shown in Fig. 4.10a. We 
see that all considered effective engineering constants depend linearly 
on the volume fraction. Clearly, vrr exhibits small oscillations in the 
considered volume fraction range. This is caused by the definition of the 
Poisson’s ratio in terms of a quotient of fitted parameters. 

To quantify the influence of changing the volume fraction, we consider 
volume fraction cg = 50% as our reference and compare the relative 
deviations of the effective elastic parameters with respect to the results 
computed for the reference volume fraction. The computed deviations 
are shown in Fig. 4.10b, where the top horizontal axis shows the relative 
deviation of the actual volume fraction cg from the reference volume 
fraction. On the vertical axis we see the corresponding relative deviations 
of the effective transversely isotropic engineering constants from the 
ones we get for cg = 50 %. The identity lines are shown in dashed gray. 


The effective Young’s and shear moduli are rather sensitive to a change 


of volume fraction. For a deviation of +2 %, the shear moduli differ by 


slightly less than +2 %, and for a variation of volume fraction by £5 %, 


the Young’s and shear moduli are changed by about +4 %. In contrast, 
the Poisson’s ratios are less sensitive to a deviation in volume fraction. 


Indeed, for +5% change in volume fraction, Poisson’s ratio is only 


changed by +2.5 % at most. 

The volume fraction is typically considered to be the microstructural 
parameter with the highest influence on the effective elastic properties of 
fiber reinforced composites. We can confirm that by our computational 
experiments. 
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4.4.2 Influence of orientation state 


In analogy to the previous section, we study the dependence of the 
effective elastic constants on the orientation. For this investigation, 
we consider all planar orientation states that may be described by 
the second-order OT A, parameterized by the diagonal component 
Aı. As already discussed in Sec. 4.3.3, these orientations vary from 
unidirectional in x-direction via planar isotropic in the z-plane to 
unidirectional in y-direction. 

We plot the effective orthotropic engineering parameters in dependence 
of A, in Fig. 4.11a. In accordance with our observations from Sec. 4.3.3, 
some of the orthotropic parameters coincide for the mentioned special 
orientation states A; = 0.0, A; = 0.5 and A; = 1.0, leading to a trans- 
versely isotropic behavior. 
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Figure 4.11: Computed effective engineering constants depending on the orientation 


A few effective engineering constants are largely independent of the 
orientation, e. g., Ey, Gyz, Gxz, whereas others change significantly, e. g., 
Ex, Ey, Vyz, Vxy. The functional dependence of the elastic parameters 
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on A, is approximately quadratic. A first glance at Fig. 4.11a suggests 
that G,, and Gyz coincide for any orientation. However, at second glance, 
cf. Fig. 4.11b, we see that both parameters actually differ, by almost 5% 
at A; = 0. However, due to the smallness of G,, and Gy, relative to 
the other engineering constants, Fig. 4.11la conveys the wrong visual 
impression. 

To investigate the orientation influence in greater detail, we limit the 
orientation states to a range of A, € [0.3,0.7], cf. Fig. 4.12. We are 
concerned with the sensitivity of the effective elastic properties w. r. t. 
small deviations of A; = 0.5. 

The effective elastic parameters for the given range are shown in 
Fig. 4.12a. The signed relative deviations are shown in Fig. 4.12b. On 
the bottom horizontal axis we see the A; component of the second-order 
OT A. We consider a planar isotropic orientation with A; = 0.5 to be 
correct. On the top horizontal axis the relative deviation of the actual 
orientation component A; from the planar isotropic value A; = 0.5 is 
shown. On the vertical axis we see the corresponding relative deviations 
of the orthotropic effective engineering constants from the ones we have 
in the planar isotropic case A; = 0.5. The identity lines are shown in 
dashed gray. 

In the vicinity of the reference orientation, the engineering constants are 
approximately linear functions of the orientation tensor component A}. 
In accordance to our observations in Fig. 4.11, an absolute deviation of 


the orientation of about +0.01, which corresponds to a relative deviation 


of +2 %, has almost no influence on certain engineering constants, such 
as the shear moduli or E,. For the remaining engineering constants 


we observe a change of up to +2.0%. Absolute deviations of +0.02 


and +0.05 lead to changes in the engineering constants of up to 43.5 % 


and +8.5 %, respectively. 


Thus, we see that the influence of the fiber orientation on the effective 
elastic tensor is of the same order of magnitude as the influence of 
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Figure 4.12: Computed effective engineering constants depending on the orientation 
for Aı € [0.3, 0.7] 


the volume fraction. This contrasts with observations for short fiber 
reinforced composites (Schneider, 2017), where the influence of the fiber 
orientation is less pronounced. 


4.4.3 Influence of matrix properties 


In this section we investigate the influence of the (isotropic) elastic prop- 
erties Em and vy of the matrix resin on the effective elastic properties of 
the SMC composites. For this purpose, the transversely isotropic elastic 
constants of the bundles need to be recomputed, as discussed in Sec. 4.3.2, 
for each new set of matrix properties. Subsequently, the corresponding 
effective stiffness of the SMC composite unit cell is computed using the 
resin moduli and the corresponding updated bundle properties. 

The effective elastic parameters resulting from a variation of the matrix’ 
Young’s modulus are shown in Fig. 4.13a. For a parameter range 
of Eu € [2.4 GPa, 4.4 GPa], the dependence of the effective Young’s and 
shear moduli on the matrix’ Young’s modulus is approximately linear, 
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cf. Fig. 4.13a, whereas the effective Poisson’s ratios depends on the 
matrix’ Young’s modulus in a nonlinear fashion. However, the relative 
deviations of the effective Poisson’s ratios are comparatively small. 


Quantitatively, an absolute change of the matrix’ Young’s modulus 


by +0.1 GPa, which corresponds to +3 %, leads to a variation of the 
elastic parameters of 0.2 % - 2.7 %. For a change by +0.2 GPa, a deviation 
of the elastic parameters of 0.4% - 5.5% is induced. For an alteration 


of +0.3 GPa in Em, which is +9 %, a deviation of the elastic parameters 
of 0.6 % - 8.5 % follows. 
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Figure 4.13: Computed effective engineering constants depending on the matrix’ isotropic 
elastic parameters 


The connection between the effective elastic parameters and the matrix’ 
Poisson’s ratio are plotted in Fig. 4.13b. Apparently, the effective 
shear moduli are almost independent of the matrix’ Poisson’s ratio, 
cf. Fig. 4.13b. In contrast, the longitudinal Young’s modulus EL 
and the Poisson’s ratios are very sensitive to a change of the matrix 
Poisson’s ratio vy, in particular close to incompressible matrix behavior 
at vm = 0.5. 
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We see that the in-plane Young’s modulus Er varies by about 0.5 GPa 
by changing vm from 0.3 to 0.4. In contrast, the out-of-plane Young’s 
modulus F; changes by 1 GPa. Thus, we see that if the in-plane Young’s 
modulus is of primary concern, 144 does not need to be determined to 
high precision. 


4.4.4 Influence of fiber properties 


In this section, we study the dependence of the effective elastic 
parameters on the isotropic elastic moduli of the fiber. Similarly to 
Sec. 4.4.3, we recompute the bundle properties when changing fiber 
properties. The sensitivity of the effective parameters w.r.t. Young’s 
modulus Er is shown in Fig. 4.14a. The effective elastic parameters 
depend approximately linearly on Er. Both the effective Young’s moduli 
and shear moduli increase for increasing fiber Young’s modulus Er, 
whereas the effective Poisson’s ratios exhibit a decreasing behavior, 
similar to Sec. 4.4.1. Compared to the previous sections, the effective 
elastic parameters are comparatively insensitive to variations of the 
fiber Young’s modulus Fr, in general. For instance, a deviation of the 


fibers’ Young’s modulus by +1 GPa leads to a deviation of the elastic 
parameters of 0.1% - 0.7%. 

Furthermore, we investigate the dependence of the effective elastic 
parameters on the Poisson’s ratio of the fibers vp in Fig. 4.14b. The 
Young’s and shear moduli remain almost constant within the entire 
considered range of vp. The effective Poisson’s ratios increase linearly 
with increasing fiber Poisson’s ratio vr. To quantify the latter observation, 


deviating from the reference Poisson’s ratio vr = 0.22 by +0.04, i.e. 
+18.2 %, leads to a deviation of the elastic parameters of 0.2% - 1.5%. 
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Figure 4.14: Computed effective engineering constants depending on the fibers’ isotropic 
elastic parameters 


4.4.5 Comparison to mean-field estimates and full-field 
computations on a pCT image 


In the first part of this section, we compare our computational multiscale 
approach to a simple mean-field estimate, as used for instance in 
Schemmann et al. (2018b) for the modeling of SMC composites. More 
precisely, we compute the effective elastic parameters based on the 
Mori-Tanaka mean-field method (Mori and Tanaka, 1973). Following the 
approach of Benveniste (1987) as presented in Brylka (2017) we use the 


expression 
ign, Fl = 
C = Cm + c& (a ( (Fo 6-0") ) + CR (Cp —Cm)* 
F 
(4.3) 
for the effective Mori-Tanaka stiffness tensor C. Here, Po = i Oe 


denotes the the polarization tensor (Ponte Castañeda and Suquet, 1998), 
čo is Eshelby’s tensor (Eshelby, 1957; 1959). The orientation 
averaging (-), can be expressed solely in terms of OTs of second- and 


where 
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fourth-order, cf. Advani and Tucker (1987). For given second-order 
OTs, the exact closure approximation, as presented in Sec. 3.4, and the 
material properties in Tab. 4.1, we compute the Mori-Tanaka estimates 
for a variation of the orientation state in analogy to Sec. 4.4.2. 


Taking a look at the effective orthotropic elastic parameters, cf. Fig. 4.15a, 
reveals a good qualitative agreement of the Mori-Tanaka predictions 
with the full-field results of Fig. 4.12a. 
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Figure 4.15: Computed effective orthotropic engineering constants based on Mori-Tanaka 
mean-field approximation for varying orientation 


Relative errors of the effective Mori-Tanaka moduli, normalized with 
the full-field homogenization results are shown in Fig. 4.15b. For the 
majority of orientations and moduli, the relative deviations do not 


exceed +7%. However, the in-plane shear modulus and two of the 


six Poisson’s ratios lead to higher errors up to +14%. In particular close 
to the planar isotropic fiber orientation state the in-plane shear modulus 
differs strongly. 
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We wish to highlight the differences in the approaches used. For the full- 
field computation, transversely isotropic bundles are embedded into an 
isotropic matrix, whereas the Mori-Tanaka approximation accounts for 
the E-glass fiber reinforcement on the single fiber level, and is insensitive 
to the formation of bundles. 

Similar to the full-field studies concerning variations of volume fraction 
as well as the matrix and the fiber properties, a comparison of the 
Mori-Tanaka mean-field results to corresponding full-field results were 
conducted. However, the results are similar in trend and magnitude to 
the study presented. For the sake of brevity, we chose not to include 
those here. 

To validate whether the bundle deposition methods described in Sec. 4.2 
captures the essential geometric features of SMC composites, we wish to 
compare the SMC unit cells generated for this work to micro-computed 
tomography (CT) scans of SMC composites in terms of the effective 
computed properties. However, the microstructure of SMC is very 
challenging in that regard. On the one hand it is necessary to consider 
a sufficiently large volume to encompass the very long fibers. On 
the other hand, for the computational procedures to work, it is not 
only necessary to separate individual fibers but to leave enough space 
between the fibers to accurately predict the mechanical deformation of 
the resin in between. Unfortunately, satisfying these two prerequisites 
simultaneously is extremely challenging. 

For the work at hand, a 1004 x 1024 x 1016 X-ray microscopy scan of 
a specimen of E-glass fiber reinforced UPPH, as detailed in Sec. 4.2.1, 
serves as the basis of our discussion. The original grey-value image 
was cropped and binarized, resulting in a 900 x 700 x 500 voxel image. 
Unfortunately, computing directly on such a binarized image is not 
possible. Indeed, a closer look at Fig. 4.1a shows that within the fiber 
bundles, distinguishing individual fibers is not possible. As a work- 
around, we approximate the bundle structure by thickening the fibers 
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of the original CT image by mathematical morphological operations 
(Haralick et al., 1987), which leads to a clustering of the fibers to bundles. 
The resulting unit cell is shown in Fig. 4.16a. The fiber thickening has 
the positive side-effect that it is possible to coarsen the resolution by a 
factor of 2, and still get almost identical effective elastic properties. Thus, 
for the succeeding, we work with a resolution of 450 x 350 x 250 voxels. 


(a) Geometry used based on a pCT scan of a SMC (b) Frobenian norm of stress (in MPa) distribution 
composite unit cell for 0.5 % uniaxial extension in x-direction 


Figure 4.16: Computation of the effective stiffness based on a pCT scan of a SMC composite 
unit cell of 450 x 350 x 250 voxels 


We analyzed the fiber orientation in the volume using the star length 
algorithm (Smit et al., 1998) applicable to binarized voxel images to 
arrive at a second-order OT of 


0.533 —0.074 —0.001 
A= | -0.074 0.427 0.008 |. (4.4) 
—0.001 0.008 0.040 


We notice that the OT is strongly diagonally dominant, indicating that the 
cell axes are well-aligned with the principal axes of the OT. Furthermore, 
the orientation is close to planar in the x-y-plane, as the zz-component 
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of A is only 4%, in agreement with visual impression. Last but not 
least, we notice that about 25 % more fibers are directed into the machine 
direction x as compared to the cross direction y. 

The thickening of the fibers leads to an artificial increase of the fiber 
volume fraction cr, which we consider as a bundle volume fraction 
cp = 64.61 %, which is larger than the cg = 50 % we use for the generated 
cells. To maintain an overall fiber volume fraction of cp = 25%, we 
adjusted the fraction of fibers within a bundle to cpg = 38.69 %, leading 
to the transversely isotropic elastic parameters for the bundles in Tab. 4.3. 


Table 4.3: Mesoscopic elastic parameters used for the uCT computation 


ERL Epr VBIT VBLT GET GBIT 
27.3 GPa 7.32GPa 0.51 0.31 2.42 GPa 2.58 GPa 


Furthermore, we use the structure tensor of Krause et al. (2010) for 
determining the direction of the bundles. The elastic parameters for the 
bundles in Tab. 4.3 and the resin properties in Tab. 4.1 serve as the basis 
for computing the effective elastic parameters of the SMC composite 
unit cell shown in Fig. 4.16a. The stress field for uniaxial loading in x- 
direction is shown in Fig. 4.16b. Bundles which are aligned in loading 
direction bear the highest stresses, whereas bundles perpendicular to 
the loading direction bear the lowest stresses. The resulting effective 
orthotropic engineering constants are collected in Tab. 4.4. 

We see that an orthotropic approximation of the effective elastic tensor 
is quite accurate. Furthermore, the diagonal components of the OT (4.4) 
translate into the relative magnitudes of the directional Young’s moduli. 
The Young’s modulus Ex in machining direction x slightly exceeds the 
transverse Young’s modulus E,. The Young’s modulus in out-of-plane 
direction E, takes a value that is only about 2/3 of Ex. The out-of-plane 
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shear moduli are almost equal, with a value about 2/3 of the in-plane 
shear modulus Gy. The values of the Poisson’s ratios lie within a range 
of 0.24 and 0.4, showing quite strong variation. 

For comparison, we generated an SMC composite unit cell and computed 
the effective elastic tensor with the multiscale approach proposed in this 
article. As we can only prescribe two-dimensional orientation states in 
the bundle deposition algorithm, we chose a diagonal two-dimensional 
orientation tensor diag(Aı, A2), s.t. the quotient A;/A2 equals the 
quotient Axx/Ayy of the first two diagonal elements of the OT (4.4), i. e., 
A; = 0.552 and Aa = 0.448. For the computations, we chose a bundle 
volume fraction of 50 % and the bundle elastic moduli of Tab. 4.1, s. t. we 
arrive at a fiber volume fraction of 25% in the end. The generated unit 
cell and an exemplary stress field according to the CT scan presented 
in Fig. 4.16 are shown in Fig. 4.17. 


(a) Generated SMC composite unit cell based on (b) Frobenian norm of stress (in MPa) distribution 
the relevant parameters of the discussed uCT scan for 0.5 % uniaxial extension in x-direction 
Fig. 4.17a 


Figure 4.17: Computation of the effective stiffness based on a generated unit cell of 
750 x 750 x 21 voxels according to the given nCT scan Fig. 4.17a 


The resulting orthotropic moduli are listed in Tab. 4.4. All six Poisson’s 
ratios and the shear moduli match almost perfectly with the full-field 
simulation results run on the pCT image. The Young’s moduli in all 
directions are slightly overestimated. 
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Furthermore, we also included the Mori-Tanaka estimates using 25 % 
fiber volume fraction and the original orientation tensor (4.4). The 
mean-field results also overestimate Ex, but recover E, quite accurately. 
However, the transverse Young’s modulus Fy is strongly overestimated. 
Similarly, all shear moduli are overestimated compared to the pCT 
results. However, the predictions of the Poisson’s ratios are quite good. 
For comparison, we included measurements of the in-plane Young’s 
moduli, obtained by two different groups. Trauth et al. (2017a) 
investigated E-glass fiber renforced UPPH resin SMC specimens, with 
a mold coverage of the produced SMC plaques of 100%. Trauth et al. 
(2017a) varied the fiber volume fraction cr between 17% — 31%. They 
observed an increase of the resulting tensile Young’s modulus E, from 
7.3 GPa to 12.5 GPa. This observation agrees with our investigations, cf. 
Sec. 4.4.1. For the comparison at hand, Trauth et al. (2017a) investigated 


specimens with a fiber volume fractions of cp = 25% + 2%. 
Furthermore, we have included the experimental investigations of 
Kehrer et al. (2018) in Tab. 4.4. They used dynamical mechanical analysis 
to determine the Young’s modulus of E-glass fiber reinforced SMC 
composite specimen. The fiber volume fraction was roughly cr ~ 25 % 
with an unknown uncertainty. 

Some care has to be taken when comparing the experimental data to 
the computational results, because the experiments are conducted on an 
entire specimen. In particular, some effects concerning the variation of 
volume fraction and the orientation are averaged out. 

The Young’s moduli E, in z-direction measured by the two groups 
almost coincide, whereas the the Young’s moduli Ey in y-direction differ 
significantly (roughly 1 GPa). However, the variances for Ey are quite 
large, so that the confidence intervals of both measurements have non- 
trivial intersection. 

All three computational approaches, full-field simulation on a uCT scan, 
the use of a generated unit cell and the Mori-Tanaka mean-field approach, 
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Table 4.4: Computed effective elastic parameters for orthortopic stiffness approximation 
and experimentally measured parameters 


pCT gener. Mori- Trauth Kehrer 
scan unitcell Tanaka (20174) (2018) 
E,inGPa | 9.42 9.55 10.02 10.96 + 0.3 10.92 + 0.6 
E,inGPa | 8.21 8.54 8.82 925410 828+05 
E,inGPa | 6.19 6.43 6.18 
Gyzin GPa | 1.95 1.95 2.03 
Gy, in GPa | 1.96 1.96 2.06 
Gyy in GPa | 3.11 3.21 3.5 
Vyz 0.398 0.388 0.407 
yz 0.368 0.365 0.379 
xy 0.342 0.341 0.341 
Vey 0.301 0.292 0.289 
Va 0.242 0.245 0.234 
Vyx 0.298 0.305 0.300 
rth +1. 0) 
E PX in% | 3.27 0.32 3.52 
cr in % 25 25 25 2542 25 


underestimate the measured Ex. Young’s modulus E, measured by 
Trauth et al. (2017a) is higher than our computed results, whereas Ey 
measured by Kehrer et al. (2018) matches quite well. 

Due to the high variation in the measured tensile Young’s moduli, 
the significance of the sensitivity studies conducted in this work is 
highlighted. At present, it is not possible to pCT scan and compute 
on an entire specimen that is used for experiments. Thus, we are com- 
paring experiments to simulation results with different microstructural 
parameters, like fiber orientation, and different length scales. Still, we 
see that for all three computational approaches, the agreement with 


experiments is reasonable. 
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4.5 Conclusions 


This chapter was devoted to extending the sequential bundle deposition 
method of Chen et al. (2018b) to be fast, accurate as well as predictive, 
and to apply the developed methods to a state-of-the art glass fiber 
reinforced UPPH based SMC composite. Understanding the effective 
behavior of SMC composites, as well as the underlying influences and 
corresponding sensitivities, is of central importance for using SMC 
composites as structural components. 

One of our key innovations concerns the use of modern orientation 
tensor closure approximations, more precisely the exact closure ap- 
proximation of Montgomery-Smith et al. (2011), and the reliance upon 
quasi-random numbers for sampling the orientation distribution. We 
found a simple formula for generating correct two-dimensional angular 
central Gaussian distribution, based on an explicit expression for the 
exact closure approximation in two dimensions. For the general, 
three-dimensional case, no such formula is known (and not expected 
to exist). The quasi-random sampling strategy enables generating high 
fidelity unit cells, where prescribed orientation data is matched closely. 
Defining the orientation tensor layer by layer permits studying the 
effects of manufacturing and placing the unidirectional bundle tapes, 
individually. 

In our opinion, the principal disadvantage of the sequential bundle 
deposition method is the disconnection of fiber bundles upon displacing 
these bundles to the next layer in case of overlap. For the elastic 
computations, this procedure did not have any negative effect. For 
inelastic computations, including fiber failure, that may no longer be the 
case. However, such an investigation is beyond the scope of this work. 
Furthermore, the curvature of the bundles (Pinter et al., 2016; Schöttl 
et al., 2021b) is not taken into account, which might have an effect. 
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In order to perform detailed sensitivity analyses on parameter influences, 
thousands of accurate microstructures with unambiguous and known 
parameter values are crucial. Therefore, these analyses cannot be 
performed on CT scans of real specimen or unit cells generated 
based on a slow or inaccurate algorithm. Based upon our improved 
SMC composite unit cell generator, we were able to conduct detailed 
sensitivity analyses quantifying the influence of all involved mechanical 
and microstructural parameters on the effective elastic properties of 
structural SMC composites. Tab. 4.5 summarizes these influences. 


Table 4.5: Sensitivity of effective elastic parameters w. r. t. input parameters. The degree of 
darkness is proportional to the sensitivity. Numbers indicate the linear correlation in %. 


As an overview, a table containing the degree of sensitivity of the 
effective elastic properties on the input parameters is included, cf. 
Tab. 4.5. For an observable quantity q, the sensitivity s w.r.t. an input 
parameter p is defined via 


q(p + Ap) — q(p) 
ap) 


q(p+ Ap) ~ q(p)[1+ sAp], i.e. = sAp. (4.5) 


The numbers shown in Tab. 4.5 are computed from the figures in this 
chapter at Ap/p = 5%, e. g., Fig. 4.10b. The coloring of the cell is directly 
linked to the correlation. The darker the cell, the higher the sensitivity of 
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the corresponding effective stiffness parameter to a change in the input 
parameter. 

We see that the elastic moduli of the matrix have a very strong influence 
on all considered elastic parameters. Also, as expected the fiber volume 
fraction is of principal importance. Similarly, the fiber orientation plays a 
major role for some selected elastic properties. Interestingly, for instance, 
the degree of fiber orientation has little influence on EL, but a strong 
influence on Er. This might explain the differences between the two 
experimental results of Tab. 4.4. 

Depending on the structural application and the subjected load, different 
elastic parameters are of principal interest. With the conducted 
sensitivity analyses and the overview Tab. 4.5 at hand, we have a 
matrix weighting influences of the microstructural and geometrical 
parameters at our disposal. Thus, we may focus on determining 
important parameters precisely, relying upon rough estimate for the less 
important ones. Detailed knowledge of these sensitivities summarized 
in Tab. 4.5 enables accounting for uncertainty in material parameters 
and measurement errors. Material science research focus, underlying 
experimental setups and measurements can be adjusted. If, e. g., an SMC 
composite part needs to elastically sustain out-of-plane loads, a reliable 
estimation of the longitudinal Young’s modulus Er is crucial. Then, 
we need to determine the bundle volume fraction as well as the matrix 
properties accurately. The fiber properties and orientation state are less 
important in this case. 

This work serves as the starting point for studying the inelastic behavior 
of SMC using more sophisticated material models, accounting for the 
thermomechanical material behavior involving stiffness degradation 
due to damage processes on the microscale. 
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A convex anisotropic damage 
model based on the compliance 
tensor 


5.1 Introduction 


5.1.1 Research contributions 


This chapter is devoted to anisotropic continuum-damage mechanics 
in the quasi-static, isothermal, small-strain setting. We propose a 
framework for anisotropic damage evolution based on the compliance 
tensor as primary damage variable, in the context of generalized 
standard models for dissipative solids. 

Based on the observation that the Hookean strain energy density of linear 
elasticity is jointly convex in the strain and the compliance tensor, we 
design thermodynamically consistent anisotropic damage models that 
satisfy Wulfinghoff’s damage-growth criterion and feature a convex free 
energy. The latter property permits obtaining mesh-independent results 
on component scale without the necessity of introducing gradients of the 
damage field. We introduce the concepts of stress-extraction tensors and 
damage-hardening functions, implicitly describing a rigorous damage- 


1 This chapter is based on the publication “A convex anisotropic damage model based on 
the compliance tensor” (Gorthofer et al., 2022b) 
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analogue of yield surfaces in elastoplasticity. These damage surfaces 
may be combined in a modular fashion and give rise to complex damage- 
degradation behavior. 

We discuss how to efficiently integrate Biot’s equation implicitly, and 
show how to design specific stress-extraction tensors and damage- 
hardening functions based on Puck’s anisotropic failure criteria. 

Last but not least we demonstrate the versatility of our proposed model 
and the efficiency of the integration procedure for a variety of examples 
of interest. 


5.1.2 Chapter structure 


We contribute to phenomenological continuum damage-mechanics with 
a tensorial damage variable. We advocate using the full compliance 
tensor as a rather natural and observable damage variable, liberating the 
engineer of the burden of selecting the appropriate damage variable 
in advance, permitting her to focus the attention on appropriate 
kinetic laws. Thus, when it comes to continuum damage-mechanics 
of phenomenological type, the proposed framework is as ab-initio as 
possible, since only the evolution of the damage surface in stress space 
needs to be identified. 

The compliance tensor has been used as the primary damage variable 
before (Ladevéze, 1983; 2002; Ladeveze et al., 2014; Baranger, 2018). Yet, 
this approach has not yet entered the main stream of damage-modeling 
frameworks. Our theoretical contributions to compliance-based damage 
models are actually twofold. For a start, we point out that the standard 
Hookean strain energy density, regarded as a function of the strain 
tensor and the full compliance tensor, is de facto jointly convex in 
both arguments. This result is surprising, and we are not aware of 
an account in the literature (although we sincerely believe that others 
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have presumably noticed this fact before without stating it explicitly, see 
Thomas and Mielke (2010) for a special case). 


Based on the compliance tensor, we develop a simple, modular 
framework for anisotropic damage mechanics. The framework provides 
the working engineer with a number of options which we believe to be of 
advantage. Indeed, due to the convexity property of the Hookean elastic 
energy, it is possible to develop a purely hardening damage-mechanics 
modeling-framework, where localization does not become an issue. 
Very much, there are materials which show a purely damage-hardening 
material response prior to sudden and brutal failure, e. g., sheet molding 
compound (SMC) composites (Fitoussi et al., 1996; 1998; Anagnostou 
et al., 2018) comprising an unsaturated polyester-polyurethane hybrid 
(UPPH) resin (Trauth et al., 2017a; Kehrer et al., 2018; Bücheler et al., 2017) 
reinforced by glass fibers (Meraghni and Benzeggagh, 1995; Schemmann 
et al., 2018a; Görthofer et al., 2019b). 

Of course, the modeling framework is not restricted to damage- 
hardening, but may be adapted to softening in a straightforward manner. 
However, the latter scenario is rather classical in continuum damage- 
mechanics, and we decided to work out the details of a hardening 
framework in the chapter at hand, essentially due to our desire to model 
SMC materials. 

To highlight the simplicity of our proposed compliance-type damage 
modeling framework, we present a first-principles development in the 
context of generalized standard models (GSMs) for dissipative solids 
(Halphen and Nguyen, 1975) and discuss the efficient resolution of the 
evolution equations in a predictor-corrector framework, see Sec. 5.2. 
Our second contribution concerns a design methodology for the 
damage surfaces which draws upon similar approaches in (associative) 
elastoplasticity (Chaboche, 2008; McDowell, 2008; Bertram, 2011), but 
takes failure criteria and multiple damage surfaces (Jin and Arson, 2018; 
Bakhshan et al., 2018; Khayyam Rayeni et al., 2020) into account. More 
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precisely, building upon Puck’s anisotropic failure criteria developed for 
continuously reinforced polymers (Puck and Schürmann, 2002; Knops, 
2008), we design specific extraction tensors and damage-activation 
functions which present a flexible arsenal of tools that, taking the 
individual damaging mechanisms into consideration, permit building 
up an accurate and fully anisotropic continuum damage model. The 
details comprise Sec. 5.3. 

For anisotropic damage models not to be judged as purely academic, it 
is of utmost importance to establish links to experimental data and 
to compare it to (dis)similar modeling approaches. In Sec. 5.4, after 
conducting computational investigations which clarify the influence of 
the different model parameters on the damage evolution and expose 
the developing elastic anisotropy upon loading, we study a plain-weave 
mesostructure of a woven carbon-fiber reinforced thermoset investigated 
by Simon et al. (2017). We show that the convex modeling framework 
permits reproducing the effective mechanical behavior of the individual 
tows and the composite quantitatively within the loading range of 
interest, see Sec. 5.4.5. 


5.2 A compliance-based anisotropic damage 
model 


5.2.1 A convex standard model for anisotropic damage 


We will describe the damage model, in a small-strain and isothermal 
setting, as a generalized standard model (GSM) (Halphen and Nguyen, 
1975), whose framework we briefly recall. In addition to the symmetric 
dx d infinitesimal strain tensor € € Sym(d), where d = 2,3 denotes the 
dimension of the ambient space, a (Banach) space Z of internal variables 
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is postulated. Furthermore, a free energy (density) 
w:Sym(d)x Z > R, (e,z)+ w(e,z), (5.1) 


a continuously differentiable function of the strain tensor e and the 
internal variables z € Z, and a force potential ®* : Z’ — [0, co], a lower 
semicontinuous, non-negative and convex function on the continuous 
dual space Z’ satisfying ®*(0) = 0, are introduced. To ensure thermody- 
namic consistency, the Clausius-Duhem inequality (CDI), see Chapter 
13 in Haupt (2002), 


where 2 = dz/dt denotes the material time derivative of the internal 
variables, needs to be satisfied. Associated to a current equilib- 
rium state (e,z) of a hyperelastic material, the Cauchy stress ten- 
sor o € Sym(d) is defined (Halphen and Nguyen, 1975; Miehe, 2002) by 


o= Ê Cez); (5.3) 


For a prescribed loading path e : [0, T] > Sym(d) on a given interval of 
time and the initial condition z(0) = zo for some zo € Z, the evolution 
of the internal variables is governed by Biot’s (dual) equation 


. 2 Ow 
ze 08 (ie 2) ‘ (5.4) 
where 0®* denotes the subdifferential of the convex function ®* 
d*(£) = {z € Z| * (E — O*(6) > (E-£)-z forall Ee 2'}, 
(5.5) 


a subset of Z’, see Borwein and Lewis (2006) for details. Due to these 
definitions, see Halphen and Nguyen (1975), generalized standard 
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materials are automatically thermodynamically consistent. Indeed, by 
Biot’s (dual) equation (5.4), 


0*(0)-O(9)> (0-9-2 for E=- le) 60 


holds. Using ®*(0) = 0, rearranging the latter inequality yields 


-ÊH (e,2).3=6.32 OE) 20, 67) 
Oz 
i. e., the Clausius-Duhem inequality (5.2) holds in view of the definition 
of stress (5.3). As it drives the evolution of the internal variables, the 
quantity € = — un (e, z) is called driving force. 


As internal variables z of our proposed continuum damage-mechanics 
model, we consider an elastic compliance tensor 


S € Sa = {S € Sym(Sym(d))|7-S[t] >0 forall r € Sym(d)\{O}}, 

(5.8) 

and a general, variable q € Q which describes the shape and size of the 
damage surface, s. t. 

z = (S,q) € Sa x Q. (5.9) 


Notice that the set S4 of (positive definite) compliance tensors is not a 
linear space. Instead, it is an open, convex subset of the linear space of 
fourth-order tensors Sym(Sym(d)) with minor and major symmetries. 
For a GSM, the CDI (5.2) will always be satisfied. However, we need to 
ensure that the (dual) Biot’s equation (5.4) guarantees that S remains an 
element of Sy, i. e., that the compliance tensor S remains positive definite. 
In contrast, the damage-surface variables we consider live in a linear 
space Q (which we deliberately keep abstract). For the specific models 
presented in Sec. 5.3, q is just a finite collection of scalar values. However, 
our arguments cover the more general case, accounting for vector- or 
tensor-valued damage-surface variables in a natural way. 
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The free energy (density) we consider is defined by 
w:Sym(d) x Sax Q>R, (e,S,q) > wele, S) + h(gq), (5.10) 


involving the Hookean elastic energy (density) 
1 
we : Sym(d) x Sa > (0,00), (e,S) => zE S!fe], (5.11) 


and an energy (density) related to the progressive degradation of the 
material, 
h:Q>R, q> hlq), (5.12) 


which we assume to be convex and continuously differentiable. 

Notice that the Hookean elastic energy we (5.11) is jointly convex in 
both variables and infinitely often differentiable. The latter property is 
immediate', as we depends on e quadratically and the Neumann-series 


representation 
(S+) =) (- "sl SeSy, LeSym(Sym(d)), (5.13) 
k=0 


valid for sufficiently small L, shows that w. is even analytic. For the 
convexity, recall that a twice differentiable function defined on an open 
convex set is convex if and only if its Hessian is positive semidefinite 
everwhere, see Theorem 3.1.11 in Borwein and Lewis (2006). A direct 
computation shows that the Hessian admits the representation 


D’w,(e,8)[£,L] = l (€ -LS [e]) -S71 [E - LS! [e], (614 


for (e, S) € Sym(d) x Sa and (€,L) € Sym(d) x Sym(Sym(d)), see Ap- 
pendix C. Any S € Sq is positive definite, and thus, the Hessian in 


1 Alternatively, representing S in matrix form, the inverse may also be represented as the 
adjugate matrix divided by the determinant, i. e., S > ST! is a rational function of the 
matrix entries. 
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equation (5.14) is non-negative. Consequently, the elastic energy is 
convex (but not strictly convex). As we assumed the energy h (5.12) to be 
continuously differentiable and convex, the smoothness and convexity 
properties of w. imply that the free energy w (5.10) is continuously 
differentiable and convex, as well. Furthermore, the formula for the 
Cauchy stress (5.3) becomes 


o = —(e,8,q) = S™' [e], (5.15) 


i. e., for a fixed compliance tensor S, the stress-strain relationship reduces 

to Hooke’s law. 

To conclude this paragraph, several remarks are in order. 

1. Using the framework of generalized standard materials for phe- 
nomenological modeling of damage is classical. For instance, Hansen 
and Schreyer (1994) study a general tensor-valued damage variable 
coupled to plasticity in such a framework, apparently unaware of the 
connection. 

2. In phenomenological continuum damage-mechanics, choosing the 
damage variable typically comes first, and the damage kinetics needs 
to be set up based on the resulting driving forces. Our approach frees 
the reader of an a priori selection of damage variable, and permits 
her to focus on the kinetics in terms of the quite natural stress-based 
driving force. 

3. The compliance tensor has the attractive characteristic that it is a 
physical quantity which can be determined experimentally. Of course, 
determining all 21 independent parameters of a stiffness tensor in 
three spatial dimensions is a daunting task from an experimental 
perspective. Still, observability of the damage variable is not ensured 
for purely phenomenological damage vectors and tensors. 

4. The compliance tensor has been used as a damage variable before 
(Ladevéze, 1983; 2002). However, its use seemed restricted to specific 
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situations, e.g., damage modeling in ceramic-matrix composites 
(Ladevéze et al., 2014; Baranger, 2018). In this work, we advocate 
using the compliance tensor as the damage variable of choice in 
greater generality. 

. It is more than well-known that the Hookean energy (5.11) is convex 
in the strain tensor. It appears much less known that the Hookean 
energy is jointly convex in the strain and the compliance tensor. When 
coupled to an energy h which makes the condensed incremental 
potential strictly convex, the resulting framework produces an 
anisotropic damage model which does not permit localization. In 
particular, associated finite-element computations are not affected by 
mesh sensitivity induced by softening behavior. We do not want 
to argue against damage localization. Rather, we wish to add a 
powerful weapon to the arsenal of continuum damage-mechanics 
when it comes to modeling stable anisotropic damage phenomena. 

. In classical small-strain elasto(visco)plasticity the (visco)plastic 
strain €, € Sym(d) serves as an internal variable. The corresponding 
stored energy (density) 


(e-&,):Cle-e,] (5.16) 


NI = 


(€, €p) > 


with a fixed stiffness tensor C = S~! is smooth and jointly convex in 
both arguments, but not strictly convex. The Hookean elastic stored 
energy function we (5.11) may be considered as a damage-analog 
of the elastic stored energy in classical elasto(visco)plasticity (5.16). 
The combined energy taking into account damage (5.11) and 
elasto(visco)plasticity (5.16) is jointly convex in all variables. If 
plasticity is neglected (e, = 0) we recover the damage case and 
for a constant stiffness (C = S7! = const.) we recover classical 
elasto(visco)plasticity. Such a model differs from the classical 
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presentation, which is typically based on either strain or energy 
equivalence (Hansen and Schreyer, 1994, Sec. 3.2.1 & 3.2.2). 


. If we regard the Hookean elastic stored energy function 


(e,C) ++ $e -C [e] as a function of the stiffness tensor C, it will not be 
convex. Indeed, its Hessian at (e, C) computes as 


(EL) + 2€- L [e] +£- C [8], (5.17) 


which may become negative (take, for instance & = e and L = -C). 
This lack of convexity is the reason why it is so difficult to design 
convex damage models for stable damage processes. Using the 
compliance tensor eradicates these issues with the help of a nonlinear 
transformation. 


. For the thermodynamics considerations at the beginning of this 


section to be valid, the “interfacial” energy (5.12) need not be convex, 
see, for example, Govindjee et al. (1995). In particular, softening 
behavior can be modeled in the compliance-tensor framework, as 
well. In that case, for obtaining a well-defined boundary-value 
problem, damage localization has to be overcome, for instance by 
adding gradient terms of the variable q to the energy (5.12). 


. The presented model cannot distinguish tensile and compressive 


loading. Indeed, the driving force T (5.18) for the compliance 
evolution computes as 
Ow(e,S, q) Ow.(e,S) 


1 
T= == = se @oa €Sym(Sym(d)), (5.18) 


which is insensitive to the involution ø +> —o. Thus, in order to 
extend our model to account for tension-compression asymmetry, the 
free energy w requires a modification, see Ladevéze (1983; 2002) and 
Ladevéze et al. (2014). 
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; m 
Damage region 1 


B#0 A 


Damaged, i F 
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B0 A | elastic region 
B=0 Unloading A B=0 


> p 
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_ 


Figure 5.1: Schematic of the admissible elastic region in (T,3)-space 


To finish presenting the two-potential model, a suitable force potential ®* 
needs to be provided, entering the evolution equation of the internal 
variables, see Eq. (5.4), 


(S, q) € 98* (T,ß). (5.19) 


In the quasi-static setting targeting a rate-independent damage model, 
we describe the force potential ®* in terms of M continuously differen- 
tiable and convex damage functions ¢; : Sym(Sym(d)) x V > R, i.e., 


0, o(T,6) <0 forall i=1,...,M, 


5.20 
+oo, otherwise. ( ) 


Such a force potential gives rise to a quasi-static damage evolution in 
terms of an elastic domain defined by the functions ¢;, in strict analogy 
to associated elastoplasticity at small strains, see Chapter 5 in Simo and 
Hughes (1998). A schematic of the admissible region based on the force 
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potential (5.20) with corresponding driving forces T and ( is shown 
in Fig. 5.1. 

However, some care has to be exercised, as the elastic domain is 
defined in terms of the compliance driving-force T, which takes the 
form T = 30 ® ø for the free energy w (5.10), in contrast to elastoplas- 
ticity, where the stress tensor ø (or a shifted version thereof) serves as 
the driving force. 


For the specific force potential ®* (5.20), Biot’s (dual) equation (5.4) 


becomes 
M M 
. 4. a¢i(T, 8) Sh. 8¢i(T, 8) 
a a and i=) hi go 62W 
involving the driving forces 
ja MED a g (5.22) 


Oq q 


for the evolution of the damage-surface variables q and consistency 
parameters p1,..., um which obey the Karush-Kuhn-Tucker (KKT) 


d , 


conditions 


To ensure that S remains in the set S4 of positive definite compliance 
tensors, a condition of the form 
9i(T, 8) 


— ap Z” forall i=1,...,M, (5.24) 


on the damage functions ¢; is sufficient. The latter condition was 
established by Wulfinghoff et al. (2017) as a criterion any physically 
meaningful vectorial or tensorial continuum-damage model should 
satisfy. In our context, the compliance tensor S serves as the damage 
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variable, and Wulfinghoff’s criterion becomes “S >0”,i.e. Sis positive 
semidefinite. 

To complete describing our model, we restrict the space of damage 
variables to Q = R™, i.e., one scalar damage variable per damage- 
activation function ¢;. We define the damage-activation functions for 
i = 1,..., M, tobe 


Qi : Sym(Sym(d)) x R— R, (T,6;) > 2T- B? — ni + H;ß;, (5.25) 


involving a (fourth-order, dimension-free) extraction tensor 


i € L(Sym(d)) with minor and major symmetries, a damage-activation 


threshold oo; (analogous to the yield stress in elastoplasticity), and a 
positive parameter H; with the dimensions of stress. 


In principle, the extraction tensor B; need not have the major symmetry 


for Eq. (5.29) to make sense. In this non-symmetric case, the term B? in 


Eq. (5.25) needs to be replaced by B}"B; in terms of the transpose 


of the extraction tensor B;. However, the framework (5.25) may be 


recovered by defining B; = \/B/"B;. Thus, by restricting the extraction 
tensor to have major symmetries we do not lose generality. Furthermore, 


as we consider the variable q; to be dimensionless, the associated driving 
force 3; has dimensions of stress and the parameter H; is necessary for 
dimensional reasons. 

In any case, for the damage function (5.25) the condition (5.24) to 
fulfill Wulfinghoff’s damage growth criterion, is automatically satisfied. 


Indeed, for any i = 1,..., M, we obtain, 
TB), _ oe . 
Le trl = r: B? r] = Bile] B: lr] = 


= |B; [r] ||? >0 forall re Sym(d). 
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In addition to the damage functions, we assume a hardening-type 
damage-surface potential of power-law type 


ee 

h(q) = XN — — dir i >0, 5.27 

(a) maa ‚m (5.27) 

involving a positive, dimension-free power-law exponent m; and a 

positive hardening parameter G; with dimensions of stress. Thus, 
according to (5.22), the damage-driving forces compute as 


In view of the force potential ®* (5.20) and the driving forces T (5.18) 
and £ (5.28), there is an elastic domain in the (extended) stress space, 
described by the conditions 


| ilo] ||? <06,+GiHig™, t= lea M,; (5.29) 


where no damage occurs. As defined in equations (5.21), the evolutions 
of the compliance and the damage-surface variables are governed by 


M 
S=2) mB? and Gg =iHi, i=1,...,M, (5.30) 
i=1 


in case of an active damage system at index i — otherwise, ù; = 0 holds. 
Several simplifications are in order. First, notice that the parameters G; 
and H; only enter (5.29) as the product G;H;. As we may redefine 
G; = H; = /G;H; without changing the elastic domain (5.29), we 
assume G; = H;. Secondly, we may eliminate the consistency parameter 
from the evolution of the compliance (5.30) and integrate to get 


N 


M .(t) 
S(t) =S +20 m 2, (5.31) 
i=l t 
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where So = S(0) is the initial compliance. Thirdly, in three spatial di- 
mensions d = 3, the compliance tensor S is described by 21 independent 
parameters. The latter formula (5.31) permits us to express the current 
compliance tensor S in terms of the internal variables q. Thus, with an 
eye towards an efficient implementation, we may a posteriori eliminate 
the compliance tensor S from the model. 


Summary of compliance-based convex damage model 
(primal formulation) 


Input. Initial compliance tensor So, extraction tensors B;, 


hardening moduli H; > 0, damage thresholds oo; > 0, power- 
law exponents m; >0 (¢=1,...,M). 


Evolution equations. For given strain path 
e: [0, T] > Sym(d), find damage-hardening variables 
q : [0,7] > RM and a stress path ø : [0, T] > Sym(d), s. t. 


holds, with initial conditions q(0) = 0, and where 


filo, qi) = |B: [o] ||? — of; — H?” and 


a= (9+5 4 a a) 


Furthermore, as stated above, notice that the equations (5.30) permit us 
to eliminate the parameters ù completely. Last but not least, in view 
of the elastic domain (5.29), we may work with the damage-activation 
functions f; for i = 1,..., M, 


fi: Sym(d)x R> R, (0,4) |Bile] ||? — o8, — Hai", (5.34) 
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instead of the original functions ¢; (5.25). For convenience, we epitomize 
the key aspects of the model in the summary box above. 


5.2.2 Computational predictor-corrector framework 


In this section, upon an implicit Euler discretization in time, we discuss a 
predictor-corrector solution strategy for the model introduced in the 
previous section in strict analogy to associative elastoplasticity, see 
Chapter 2 in Simo and Hughes (1998). Suppose that a number of discrete 
time steps 0 = to < tı <...<ty_1 < tyn = T is given, together with 
prescribed strain tensors eo, €),...€y, an initial compliance tensor So 
and the initial damage-hardening variable qo = 0 € RM. For any 
n=0,...,N — 1, dropping the subscript n + 1 for simplicity of notation, 


, 


we seek (o, q) € Sym(d) x R™ solving the system of equations 


M 
qi 2 
e= |S% +25 = ) [o] 
( a Hi 


filo,gi) <9, qi— qin 20, (qi— din) file, gi) =0, i=1,...,M, 
(5.35) 


with the damage functions f; (5.33);. With a computational resolution 
in mind, we rewrite the system (5.35) in terms of active sets. For any 
(o, q) € Sym(d) x RM, the active set A(o,q) is defined as 


Alo,g) = {i € {1,2,..., M} | filo, a) = 0}, (5.36) 


collecting all indices of inequality constraints that are either violated or 
satisfied exactly. Then, as a consequence of the complementarity condi- 
tion in the system (5.35), (7, q) € Sym(d) x R™ solves the system (5.35) 
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precisely if it is satisfies q; > in (i =1,...,M) and solves 


M 
qi 2 
So + 2 — B; | [ol] =e 
(s IH, J (5.37) 


filo,q) =9 forall ie Alo,g). 


We solve the latter problem by an active set strategy (Bergounioux 
et al., 1999; 2000), i.e., by solving the system (5.37) with a Newton 
method, updating the currently active set at each Newton iteration and 
accounting for the constraints qi > gin (i = 1,..., M) via backtracking. 
The details comprise Alg. 1, where y € (0,1) is a backtracking factor. 
We use a backtracking factor of y = 0.9 in our presented examples, 


see Sec. 5.4. 
ott 
Elastic 
predictor 
fa=0 region due to damage 
Admissible region at tp 
finti =0 
f Lae 0 
fz =0 


fon = fon =0 


Figure 5.2: Evolution of the elastic region upon loading within a predictor-corrector 
framework 


As long as the damage constraints are linearly independent, due to the 
established connections of active set strategies to semi-smooth Newton 
methods, see Hintermiiller et al. (2002), a locally superlinear convergence 
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behavior can be expected. A schematic of the predictor-corrector strategy 
is shown in Fig. 5.2 with 


2 
residual(o,q) = ene where 
oO 
M r = 
Rı=o- (s + > T, ) [e], and (5.38) 


M 
Rə = ) max{0, fi(o, qi)} 

i=1 
for measuring convergence. Whenever the trial stress fails to be 
contained in the elastic region, an iterative process is initiated which 
ensures that the final stress state again lies on the boundary of the elastic 
domain. For the latter, both the elastic region may grow - as a result 
of the damage-hardening - and the stress may decrease — owing to 
increasing compliance. 
For solving problem (5.37), we assemble the Newton system for the 
active set A(o,q) 


= e- (5+22%, B?) jo} 


1 , 
= -- (|B: [ø] |? - o8; - Hai") 


Hi 
(5.39) 
for alli € A(o,g), where we divided the second line by H; to ensure a 


symmetric Newton system. 
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Algorithm 1 Predictor-corrector strategy (e,g;.n) with model parameters 


(So, 


i, Hi, 00,,,m;) and algorithm parameters (maxit, tol, y) 


m 


O O0: Sa). (OS ON et a 


PrP PrP RP RP eB 
Sup Oy ES 


: Elastic predictor 


17 An 


: if all f; < 0 then 


no damage evolution, elastic predictor step correct V 
else 
Damage corrector 
kel > Iteration counter 
Update residual (5.38) 
while k < maxit and residual > tol do 
A+ A(o,q) 
assemble Newton system (5.39) 
solve for (Ao, Aq) > Aq :=0forif A 
scl > Full step size 
(o,q)- (o +s Ao,q4+sAq) 
j+0 > Counts backtracking steps 
residual,jqg < residual 
Update residual (5.38) 
while residual > residualoıa or qi < din for some i do 
(Backtracking, typically y = 0.9) 
(0,9) + (o + (ys — 8)Aa,q + (ys — s)Aq) 
sys > Reduce current step size 
Update residual (5.38) 
jogjgtl 
end while 
k¢ek+1 
end while 


: end if 

: compute Caigo 

: Output 

: return a, q, Caigo 
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5.3 Damage models with Puck-type extraction 
tensors accounting for average stresses 


5.3.1 Basic idea 


Puck introduced strength-estimation models for composites reinforced 
by continuous fibers (Puck and Schiirmann, 2002; Knops, 2008) based 
on specific failure scenarios that are commonly observed in post-critical 
investigations of failed specimens. For the current chapter at hand, we 
will use these so-called Puck cases as primary drivers of the anisotropic 
damage evolution presented in Sec. 5.2. More precisely, we will 
investigate the Puck cases individually, and determine proper extraction 


tensors (Br - Bry) in sections 5.3.2 - 5.3.5. These Puck-type extraction 


tensors are motivated by the stress and corresponding damage states 
present in the fiber bundle mesostructure of sheet molding compound 
(SMC) composites (Dumont et al., 2007; Gorthofer et al., 2019b) in an 
averaged setting. 

We introduce a local Cartesian coordinate system {e1, ea, e3}, s.t. the 
fibers are aligned to the eı-direction, see Fig. 5.3a. Then, the stress state o 
may be decomposed into blocks 


O11 | 012 013 


g= 012 | 022 023 ; (5.40) 


013 | 023 933 


where 01; is the stress in fiber direction, the lower right block describes 
the stresses in the plane orthogonal to the fiber direction, and (012, 013) 
collects the remaining shear stresses. Adopting ideas of Puck (Knops, 
2008; Puck and Schiirmann, 2002), we distinguish four basic cases which 
drive the damage evolution in a fiber bundle, for instance. 
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(I) Normal loading in fiber direction: 014 
Figs. 5.3a and 5.3b 

(II) Normal loading perpendicular to fiber 
Figs. 5.3c and 5.3d 


direction: 022, 033 


(III) Shear loading perpendicular to fiber direction: 023 


Fig. 5.3e 


(IV) Shear loading in fiber direction: 012, 013 


Fig. 5.3f 


Damaged region 


Fibers BZ 


tion 


(d) Compression _ to fiber direc- (e) Shearing | to fiber direction (f) Shearing in fiber direction 


tion 


Figure 5.3: Regions of major damage (blue) resulting from different loading scenarios in a 


cell with aligned fibers (dark green) 


The loading scenarios shown in Fig. 5.3 are only examples, e.g., 


loadings perpendicular to the fiber direction need not necessarily follow 


direction ea. Instead, any other direction in 


the e2-e3-plane could be 


used, as well. Nevertheless, we may regard a general loading scenario 


as a superposition of the four introduced cases. In the following sections, 


we will derive appropriate extraction tensors ( 


By - Brv) corresponding to 


101 


5 A convex anisotropic damage model based on the compliance tensor 


each of the four presented cases based on averaged stress conditions. The 
presented model cannot distinguish between tensile and compressive 
loading, as the driving force T (5.18) is quadratic in the stress o. 
Consequently, the six sketched loading scenarios in Fig. 5.3 reduce to the 
mentioned four cases, as the scenarios Fig. 5.3a and Fig. 5.3b, as well as 
Fig. 5.3c and Fig. 5.3d coincide for our model. 


5.3.2 Normal loading in fiber direction 


The first damage case is governed by loading in fiber or bundle 
direction, respectively, and thus solely concerns the stress 11. For fiber 


direction eı, the fourth-order extraction tensor By extracting the stress in 


bundle direction c11 from an arbitrary stress state ø is given by 


(=e: (5.41) 


The associated damage function (5.33), with case I extraction tensor By 
(5.41) reads 
fi(a,g) = oii — H’q™ — 0%, (5.42) 


and will solely induce a decrease in the Young’s modulus in e,-direction. 
As a general note, although the damage parameters like oo, H, m may 
differ for the considered cases I to IV, we do not include additional 
subscripts for the sake of readability. 


5.3.3 Normal loading perpendicular to fiber direction 


To quantify damaging due to normal loading in any direc- 
tion 9? 5 n L e; perpendicular to the fiber direction (a unit vector 
with S? = {x € R? | ||x|| = 1}), for a general stress state o, we measure 
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the average normal stress perpendicular to the fiber direction 


1 20 gi 

— 4 
an, n(6)** [o] dé (5.43) 
with n(0) = (0, cos 0, sin 0). The latter average may be represented in the 


form 
Í 2m 


n(0)®* [o] d0 = By [ø] (5.44) 


2m Jo 


in terms of the extraction tensor By; 


n(0)*4 40 (5.45) 


with n(0) = (0, cos 6, sin 0). 


We evaluate the integration for each component of the extraction tensor 


separately. As n L eı, all components of By; with at least one index "1" 


are zero. The remaining components are 


27 
BE gcse cos! 0 d8 = 2 (5.46) 
2222 Ir å 8’ . 
1 27 
Bigs 5 f cos? 0 sin? 0 d0 
1 
B3322 B3323 B3332 Byp32 Bias g’ (5.47) 
2m 
Bu sinf 6d@ = : (5.48) 
3333 Ir ð g’ : 


1 27 
Bin = 37 i cos? Ø sin 0 d0 


B3332 = B3322 = B3922 = 0, (5.49) 


1 27 
Bg =| cos 0 sin? 0 d0 


B3523 = B5hs3 = Bass = 0. (5.50) 
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Hence, the extraction tensor is 


1 1 1 
T= je + eS?) 22 + ACH = 8?) 23 5(€2 @s e3)®?. (5.51) 


This extraction tensor By; is identical to the fourth-order fiber-orientation 


tensor for a planar isotropic orientation, see Advani and Tucker (1987). 
The composition of the extraction tensor with itself is 


1 1 1 
i= glee) + eg”)? + ao (ef? — ef")? + (er @ses)®?. (5.52) 


The damage function (5.33), involving the case II extraction tensor By 
(5.51) reads 


1 
fulo, q) = 32 (5032 + 5033 + 6022033 + 4033) — H7q’™ — og. (8.53) 


5.3.4 Shear loading perpendicular to fiber direction 


In addition to damage caused by normal loading, we also want to 
account for shear-loading induced damage. For a general stress state a, 
the average shear stress transverse to the fiber direction is given by 


1 2T 
— | (a0) 8s m(6))*? [o] ad (5.54) 

2T 0 
with n(@) = (0,cos6,sin6) and m(@) = (0,—sin6,cos@). We may 
rewrite this expression 


1 27 


(n(0) 8s m(0))®? [o] dO = Bm [o] (5.55) 


2m Jo 


in terms of the extraction tensor By; as 


1 


Qn 


I T (m(b) 8s m(0))2? d0 (5.56) 
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with n(0) = (0, cos 6, sin 0) and m(0) = (0, — sin 0, cos 6). 

We evaluate the integration for each component of the extraction tensor 
separately. Asn L m L e; and nı = mı = 0, all components with at 
least one index "1" are zero. For the remaining components we get 


wo 1 
Ei, = f cos? 0 sin? 0 d0 = 7 (5.57) 


1 27 
Bis = f sin? 0 cos? 0 d0 


1 
= By = — 3: (5.58) 
Il 1 ae 4 4 2 2 
B3323 = =] = (sinf 0 + cos“ 0 — 2sin” 0 cos” 0) dé 
QT 0 4 
1 
B3532 = B3332 = B3223 8° (5.59) 
II a me r 1 
B3333 = on 5 sın 6 cos 0d0 = g’ (5.60) 
1 27 
Des / (sin? 0 cos 8 — sin 8 cos? 0) dé 
At Jo 
B3332 = B3522 = B3322 = 0, (5.61) 
1 27 
B5332 = a, (sin 0 cos? 0 — sin? 0 cos 0) d0 
B3323 = B3333 = B2533 = 0. (5.62) 
The extraction tensor for case III therefore has the form 
1, 82 3282 , 1 @2 
Bin = 8 (e5 — eş ) + „(e2 Qs e3)°”. (5.63) 
The composition of the extraction tensor with itself is 
2 1 
m = go. (5.64) 
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The damage function (5.33), for the case III extraction tensor By (5.63) 


reads 


1 
fulo, q) = 35 (032 + 033 — 2022033 + 4033) — H’q" — ob. (5.65) 


Comparing the extraction tensors By (5.51) and By (5.63), some 


similarities between these tensors become apparent. In fact, these 
similarities reflect the relationship between normal loadings (especially 
compression) and shear loadings perpendicular to the fiber direction, 
that are familiar from undergraduate engineering mechanics (Hibbeler, 
2001), i. e., Mohr’s circle (Mohr, 1900). 


5.3.5 Shear loading in fiber direction 


To evaluate damage induced by shearing in fiber direction eı, we 

evaluate 7 
1 T 
= (n(0) @s eı)” [o] dd (5.66) 
27 Jo 

with n(@) = (0, cos 6, sin 0), which we represent in the form 


1 27 


22 iø] dd = Bıv[e] (5.67) 


(n(6) @s e1) 


2T Jo 


with extraction tensor Bry 


1 27 
Biv = 5 / (n(0) 8s e1)®? d0 (5.68) 
2T Jo 


and n(0) = (0, cos 0, sin 0). 

Again, we evaluate the integration for each component of the extraction 
tensor separately. The only non-zero components are B2121, B2112, B1221, 
Bıaı3, B3131, B3113, B1331 and B1313. These components are computed 


106 


5.4 Computational investigations 


as 


Í 2r 
B i= =| cos” 0 d0 


1 
Bili = Blyn = BY = $, (5.69) 
1 2T 
B= ae f sin? 6 d0 
Bins = Bia = Biss = g (5.70) 
The resulting extraction tensor for case IV may be expressed via 
1 go, 1 22 
Iv = z(e Ss e2) + „ie Ss e3) 3 (5.71) 
The composition of the extraction tensor with itself is 
2 1 
Iv 5 v8 (5.72) 


Inserting case IV extraction tensor Bry (5.71) into damage function (5.33), 


yields 


fv(o,g) = (012 + 013) — H’g" — 95. (5.73) 


1 
8 
5.4 Computational investigations 


5.4.1 Setup 


We integrated the proposed damage model as a user-defined subroutine 
into an in-house OpenMP-parallel FFT-based computational homoge- 
nization code written in Python 3.7 with Cython extensions (Behnel et al., 
2011) and FFTW (Frigo and Johnson, 2005) bindings, as described, e. g., 
by Schneider (2018). The balance of linear momentum was discretized 
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ona staggered grid (Schneider et al., 2016b) and the ensuing nonlinear 
systems of equations were solved by a Newton-CG scheme (Gélébart and 
Mondon-Cancel, 2013; Kabel et al., 2014; Wicht et al., 2020). Furthermore, 
we implemented Alg. 1 as a C++ subroutine (Stroustrup, 1999) in 
the commercial Finite Element (FE) software ABAQUS/Standard by 
Dassault Systèmes Simulia (Dassault Systemes, 2014). The computations 
ran on 6 - 12 threads on a desktop computer with 32 GB RAM and 
an Intel i7-8700K CPU with 6 cores and a clock rate of 3.7 GHz. The 
plain-weave presented in Sec. 5.4.5 was computed on a workstation with 
two AMD EPYC 7642 and 48 physical cores each, enabled SMT, and 
1024 GB of DRAM. 


For the following studies, we use the isotropic elastic parameters of 
unsaturated polyester polyurethane hybrid (UPPH) resin and E-glass 
fibers, see Kehrer et al. (2018), respectively, if not specified otherwise. 
Furthermore, we use the damage parameters 09, H and m listed in 
Tab. 5.1 as a point of departure for parameter variation and the different 
introduced damage cases. Due to the small-strain setting, we limit the 
strain axes to 5 % in magnitude. 


Table 5.1: Standard material parameters (Kehrer et al., 2018) and reference damage 
parameters, serving as point of departure depending on the corresponding damage case 


UPPH matrix E-glass fibers Damage parameters 


Ey =34GPa Ep=72GPa ou € [5,30] MPa 
um = 0.385 Vp = 0.22 H € [30,80] MPa 


mie 1 
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5.4.2 Numerical studies on integration-point level 


Parameter study for Puck-type extraction tensor I 


The first study concerns the effects of the damage parameters oo, H 
and m on the stress and damage evolution. For this purpose, we 
investigate the model behavior for one active damage function and the 
Puck-type extraction tensor I (see Sec. 5.3.2) only. We vary the parameters 
and evaluate the stress-strain curves, as well as the normalized Young's 
modulus E11/E?, for uniaxial extension in eı-direction and a prescribed 


strain £11. 


For the case at hand, we extract the current Young’s modulus in eı- 
direction from Eq. (5.31) by 


HEN, 


igs  EE 
11 H +2qE9,’ 


(5.74) 
where EP, stands for the initial Young’s modulus in eı-direction. Please 
note that Ej, may be regarded as a component of a (fourth-order) tensor, 
as it is the 1111 component of the corresponding compliance S. For the 
sake of simplicity, we will use such shorthand index notation for certain 
Young’s moduli throughout the remainder of this work. 

The influence of the damage-activation threshold oo on the stress and 
normalized stiffness in eı-direction is shown in Fig. 5.4. The higher the 
damage-activation threshold oo, the later the damage evolution initiates 
w.r.t. the applied strain £11. In Fig. 5.4a, we observe damage to initiate as 
soon as the stress a1; equals the damage-activation threshold oo, which 
is expected. The convex hardening nature of our model gives rise to a 
decreasing slope of the stress-strain curve. This slope tends to zero at 
infinity, but remains non-negative. 

With increasing damage-activation thresholds oo, the stress-strain 
behavior approaches a plateau beyond the elastic region, in which 
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— o = 10 MPa — on = 30 MPa — oo = 50 MPa — oo = 70 MPa — oo = 90 MPa — oo = 110 MPa 
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(a) Stress-strain curve (b) Normalized Young’s modulus vs. strain 


Figure 5.4: Varying the initial stress oo for the proposed model 


an increase of strain does not induce a further increase of stress. The 
reduced stiffness £1; in eı-direction equals the slope of the stress-strain 
curve during unloading (which returns to the origin in our pure elastic- 
damage framework, see Sec. 5.4.2). Whereas the plateaus are more 
pronounced for higher damage-activation thresholds ao, the increase in 
damage and the (normalized) stiffness reduction are less pronounced, 
see Fig. 5.4b. 


Due to the thermodynamically consistent GSM framework of our 
proposed model, an upper bound for the damage variables is ensured, 
governing the asymptotic behavior of the (normalized) stiffness reduc- 
tion, see Fig. 5.4b. Evaluating the CDI (5.2) for the considered case at 
hand, we find the upper bound for the damage variable w.r.t. Puck 
case I tobe q < “oj, /H?. 

The effect of changing the hardening parameter H on the stress and 
damage evolution is shown in Fig. 5.5. As the damage-activation 
threshold co remains unchanged for this study, the elastic regions (011 
below 30 MPa) coincide for all stress-strain curves, see Fig. 5.5a. 
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— H = 10 MPa — H = 70 MPa — H = 130 MPa — H = 190 MPa — H = 250 MPa — H = 310 MPa 
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Figure 5.5: Varying the hardening parameter H for the proposed model 


In the damage region, the slope of the stress-strain curve increases 
with the hardening parameter H. Indeed, the hardening parameter H 
describes the hardening capacity of the model. For H — 0, the slope 
tends to zero and approaches the plateau already observed in Fig. 5.4a. 
For H — oo, an active damage function f (5.33); will become inactive 
for small values of the damage variable q — 0, resulting in a damage 
region that can hardly be distinguished from a purely elastic material 
behavior. 

Independent of the hardening parameter, damage evolution initiates at 
the same strain level (and therefore the same stress level), see Fig. 5.5b. 
The stiffness reduction is inversely proportional to the hardening 
parameter H. 

In Fig. 5.6, we fix the damage-activation threshold oo as well as the 
hardening parameter H and vary the power-law exponent m. In contrast 
to the influence of the hardening parameter H (see Fig. 5.5), the slope of 
the stress-strain curve is inversely proportional to the exponent m in the 
damage region, see Fig. 5.6a. 
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— m = 0.2 — m = 0.6 — m = 1.0 — m = 1.4 — m = 1.8 — m = 2.2 
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Figure 5.6: Varying the power-law exponent m for the proposed model 


For increasing exponents m, the stress-strain curves approach the 
plateau-like behavior. For small values of m, after exceeding the damage- 
activation threshold oo, the stress-strain curves remain approximately 
linear and only develop a significant amount of damage for higher 
loading levels. 

Fig. 5.6b shows the damage evolution to be inversely proportional to 
the exponent m, leading to a lower remaining (normalized) stiffness 
component F11 for higher exponents m. 

For representing the stiffness tensor C =S~! graphically, we use 
the Young’s modulus surface (YMS) plots introduced by Böhlke and 
Brüggemann (2001), i.e., for fixed compliance tensor S € S4, the image 
of the nonlinear mapping 


S? > R’, nv E(S,n)n, (5.75) 
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where the Young’s modulus E(S, n) in direction n is determined by 


1 


nz (ngn) Singen] 


(5.76) 


Asymmetry properties of the stiffness tensor C become apparent in the 
corresponding YMS plot. Examples of such YMS plots are shown in 
Fig. 5.7. 


(a) Initial isotropic state (b) Final anisotropic state 


Figure 5.7: YMS plots (see Böhlke and Brüggemann (2001)) showing the reduction of the 
stiffness tensor based on Puck case I 


The initially isotropic stiffness tensor with UPPH material parameters 
(see Tab. 5.1) has a spherical shape, as shown in Fig. 5.7a. As discussed 
above, the induced damage model based on Puck case I leads to a 
reduction of the stiffness Eıı in the eı-direction. The YMS plot is 
contracted in the direction of loading, whereas the directional Young’s 
moduli in the orthogonal plane remain unaffected, see Fig. 5.7b. 
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Stiffness reduction for different Puck-type extraction tensors 


In the following, we will discuss possible damage evolutions and 
corresponding stiffness reductions based on Puck-type extraction ten- 
sors II, HI and IV (see Sec. 5.3). For these cases, the influence of the 
damage-hardening parameters oo, H and m is similar to case I, which 
we discussed in Sec. 5.4.2. We consider an initially isotropic stiffness 
tensor with UPPH material parameters. The corresponding YMS plot is 
shown in Fig. 5.7a. Specific loadings will evoke a damage evolution due 
to the Puck cases IT, II and IV. We apply a normal strain £22 for case II, 
a shear strain £23 for case III and a shear strain £13 for case IV, forcing 
the respective complementary stress components to zero. The resulting 


YMS plots are shown in Fig. 5.8. 


a sere 
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Tieder g a a 
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(a) Case II (b) Case II (c) Case IV 


Figure 5.8: YMS plots (see Böhlke and Brüggemann (2001)) illustrating the reduction of 
the stiffness tensor based on the Puck cases II, M and IV 


For all three cases the Young’s modulus in eı-direction remains 
unchanged. As shortly discussed in Sec. 5.3.4, Puck cases II and III 
are interlinked due to similar effects of averaged normal loadings and 
shear loadings perpendicular to the fiber direction. Inspecting Figs. 5.8a 
and 5.8b, corresponding to Puck cases II and III, we observe a reduction 
of the Young’s moduli within the e2-e3-plane for both cases, but with 
different characteristics. For the loading scenarios considered here, the 
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stiffness reduction in directions e2 and e3 is more pronounced for Puck 
case II compared to Puck case III. Fig. 5.8c shows that damage based 
on Puck case IV does not affect the Young’s moduli in the e2-e3-plane. 
Young’s moduli in the e;-e3-plane and e-e2-plane are reduced equally, 
leading to a transversely isotropic stiffness with the fiber direction eı as 
the axis of symmetry. 


Non-monotonic loading 


To show the capabilities of our model in general, we perform loading- 
unloading experiments for different loading directions in a successive 
fashion. To mimic uniaxial normal loadings and corresponding shear 
loadings, we subsequently apply normal and shear strains £11, €22, €33, 
€23, €13 and £12, each with mixed boundary conditions permitting solely 
the corresponding stresses 011, 022, 033, 023, 013 and 012, to be non-zero, 
see Fig. 5.9a and Fig. 5.9b. Lateral contraction is permitted. Between 
each of these six loading steps we unload to zero strain and stress. 


Stress oj in MP 
S 


En — En — Fs — Em — EN — En Ti — 222 — Tss — 735 — 015 — 13 a au qur av 


(a) Strain vs. time (b) Stress vs. time (c) Damage evolution vs. time 


Figure 5.9: Complex loading history addressing each stress-tensor component separately 


The resulting evolution of the stress components is plotted in Fig. 5.9b. 
We see a linear elastic region for each individual loading step and a 
damage region for all but the 733 and o12 cases. On these occasions, the 
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threshold for damage initiation is not triggered. As all normal loadings in 
the e2-e3-plane evoke damage due to case II, 733 cannot induce damage 
beyond the level previously induced. The stress-strain curves for the 
shear stresses 013 and c12 exhibit a similar behavior. These observations 
are also reflected in the damage evolution, see Fig. 5.9c. Due to the 
applied stress c22, the damage variable representing case II increases. 
For increasing stress 033, the damage variables remain unaffected up to 
the level 033 = 022. The same line of argument applies to o13 and 012. 
Again, we observe the connection between Puck cases II and IL i.e., 
each of the damage variables increases whenever a loading scenario is 
applied which progresses the other case. 


(c) Loading step 3 (€11, €22 and 
€33) 


(d) Loading step 4 (€11, €22, €33 (e) Loading step 5 (€11, €22, E33, (f) Loading step 6 (€11, €22, €33, 
and £23) €23 and € 3) E23, €13 and £12) 


Figure 5.10: Evolution of an initially isotropic stiffness upon complex loading, see Fig. 5.9a, 
visualized via YMS plots 
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The damaged stiffness tensors corresponding to each loading step are 
visualized via their YMS plots in Fig. 5.10. The different colors of the 
plots represent the different loading steps shown in Fig. 5.9. Note that 
the ranges of the axes are adjusted accordingly and therefore vary from 
plot to plot as the damage increases from step to step. 

Comparing Fig. 5.10b and Fig. 5.10c, as well as Fig. 5.10e and Fig. 5.10f, 
we see that, in accordance with Fig. 5.9c, no further damage is induced 
between these loading steps. The presented YMS plots demonstrate the 
capability of our model to evolve the stiffness tensor in a complex and 
anisotropic way. The model is capable of handling any initial stiffness, 
not restricted by a specific symmetry class, i. e., transversely isotropic 
or orthotropic. Furthermore, the stiffness tensor may also develop 
anisotropy — within the permissible set S4 (5.8) — as a result of a damaging 
process, owing to the introduced damage functions. 


Multiaxial loading with increasing loading level 


In a study combining the six loading steps from the previous section 5.4.2 
to one superposed loading case, we apply a three-dimensional strain 
state s. t. all strains are active simultaneously. To investigate the model 
behavior, the predicted stiffness degradation, as well as the evolution 
of the damage variables more closely, we gradually increase the strain 
levels in five steps from 0% to 5% via 1% increments with intermediate 
unloading, see Fig. 5.11a. 

As in Sec. 5.4.2, we analyze the stress and the damage evolution, as well 
as the stiffness reduction. Fig. 5.11b shows the evolution of the individual 
stress components during the combined loading. We see that after 
each loading-unloading step the level of damage increases. This is also 
reflected in the evolution of the damage variables, see Fig. 5.11c. After 
each loading-unloading step, the damage variables continue to increase 
whenever the maximum stresses of the previous step are exceeded. 
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Stress oi in MPa 
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(a) Strain vs. time (b) Stress vs. time (c) Damage evolution vs. time 


Figure 5.11: Step wise increase of multiaxial loading to evoke all stresses and damage 
functions simultaneously 
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(a) Strain 0% (b) Strain 2% (c) Strain 5% 


Figure 5.12: Evolution of an initially isotropic stiffness during multiaxial loading steps, 
visualized via YMS plots 


The YMS plots corresponding to 2% and 5% loading, as well as the 
initially isotropic YMS plot, are shown in Fig. 5.12. The stiffness 
tensor gradually reduces based on all four introduced Puck cases (see 
Sec. 5.3) simultaneously. Apparently, the multiaxial loading evokes a 
stiffness reduction in all directions, as we observe a superposition of 
the individual loading scenarios investigated in Sec. 5.4.2 and Sec. 5.4.2. 
Please note, that, similar to the previous section, the ranges of the axes 
are adjusted accordingly and therefore vary from plot to plot. 
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Cyclic tensile loading with increasing loading level 


We conduct an analysis of the model response upon cyclic loading via 
uniaxial extension. We successively apply a normal strain £11 to induce 
cyclic uniaxial loading in the e,-direction. We apply the strain £11 in five 
cycles from 0% to 5% with an increasing magnitude of 1% per cycle. 
For this analysis, we restrict to Puck case I. The resulting stress-strain 
curves and damage-strain curves are shown in Fig. 5.13. 


— Cycle 1 — Cycle 2 — Cycle 3 — Cycle 4 — Cycle 5 
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Figure 5.13: Cyclic tensile loading with increasing loading level for Puck case I 


Upon loading, the material behaves linear elastically until a specific 
critical stress threshold or the maximum stress of the previous cycle is 
reached, see Fig. 5.13a. During the in-between unloading to £11 = 0% 
and reloading, the damage variables do not evolve further. Besides, the 
pure damaging character of the model is highlighted, as no remaining 
residual strains occur. The step wise evolution of the damage variable q 
is shown in Fig. 5.13b. 
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Note that the presented model is capable to predict damage onset upon 
both, tensile and compressive loading. Due to the definition of the 
damage functions the quadratic nature of the driving force T (5.18), 
the damage evolution is driven in similar fashions by both, tensile and 
compressive stresses. Considering a single damage variable q for, both, 
the tensile and the compressive regime, would lead to a combined 
damage evolution. Accounting for tension-compression asymmetry 
requires an extension of the model at hand, see the conclusion. 


5.4.3 Response for a continuous fiber microstructure 
Application of separate loading cases 


After discussing the model response for homogeneous stress states, 
we account for heterogeneous stress states in two ways to show the 
basic feasibilities of our damage model. First, we shall investigate 
a microstructure with a continuous fiber reinforcement. In the next 
sections, we will turn our attention to mesoscale simulations. 

To account for damage evolution in the matrix, we introduce two 
extraction tensor corresponding to the spherical and deviatoric projectors 
of fourth-order 


1 
sph = P; = z781, dev = Po = IS — P3, (5.77) 


that allow describing a damage evolution in response to dilatation 
and distortion. We use damage-activation functions based on these 
two extraction tensors (5.77) and corresponding damage parame- 
ters oo = 30 MPa, H = 130 MPa and m = 1 for both cases. Furthermore, 
the matrix of the fiber reinforced microstructure is endowed with 
the elastic properties of UPPH, as defined in Tab. 5.1. The fibers 
are modeled in a purely elastic fashion using the elastic moduli of 
E-glass, see Tab. 5.1. The continuous-fiber reinforced microstructure 
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is geometrically modeled by 50 inclusions with a diameter of 13 um and 
a total volume fraction of about 40 %. We generated the microstructure 
by the adaptive shrinking-cell algorithm of Torquato and Jiao (2010). 
The setup represents aligned fibers in a UPPH matrix as present in SMC 
composite bundles (Meyer et al., 2020; Kim et al., 2011b), where the fiber 
direction e; coincides with the primary direction of the Puck cases. 

In three different loading scenarios, we apply three different macroscopic 
strains via mixed boundary conditions, see Kabel et al. (2016) for 
details. For each scenario, we analyze the induced damage fields of 
the associated variables qsph and qaev. In scenario 1, we apply the 
macroscopic normal strain €22 perpendicular to the fiber direction (in 
horizontal direction). In scenarios 2 and 3, we apply macroscopic shear 


strains £19 and €93, in fiber direction and transverse to the fiber direction. 


sph 
0.15 
0.10 
0.05 
0.00 
(a) Loading E22, dsph (b) Loading €12, dsph (c) Loading €23, dsph 
Idev 


(d) Loading 22, qdev (e) Loading E12, qdev (f) Loading E23, dev 
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Figure 5.14: Model response for a continuous-fiber reinforced microstructure based on 
spherical and deviatoric damage and three different loading cases 
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The average runtime for a resolution of 256 x 256 pixels and 50 time 
steps was about 100 s on 12 threads. An accompanying resolution study 
is discussed in Sec. 5.4.3. Fig. 5.14 shows the damage fields for qspn and 
qaev for the introduced cases and corresponding to the different loading 
scenarios. 

Fig. 5.14a shows that normal loading perpendicular to the fiber direction 
leads to a dilatation-triggered damage evolution in the respective 
direction, as a result of stress concentrations at the inclusion boundaries. 
As a consequence of the complexity of the induced stress state, damage 
due to distortion is initiated, as well, see Fig. 5.14a. In general, the 
stress level is higher for regions with more closely packed inclusions, 
inducing a significantly higher level of damage in those regions. Damage 
initiates at the inclusion boundaries and evolves in the loading direction, 
deflected by other inclusions. Shear loading in fiber direction leads to an 
associated damage evolution due to distortion, as shown in Fig. 5.14e. 
As the applied shear is oriented in fiber direction, the deformation is 
not hindered by these fibers and spherical stresses do not occur. Hence, 
damage due to dilatation does not evolve, see Fig. 5.14b. Applying 
a macroscopic shear perpendicular to fiber direction evokes both, the 
evolution of damage due to spherical stresses, see Fig. 5.14c, as well as 
deviatoric stresses, see Fig. 5.14f. 


Resolution study 


The presented resolution study demonstrates that the proposed dam- 
age model leads to mesh-independent results even without gradient 
enhancement. This does not come by surprise, as we specifically 
designed such a hardening-type damage model. Still, even in the case of 
hardening, a resolution study is imperative to ensure mesh-independent 
results. In particular, we will justify the resolution employed in Sec. 5.4.3. 
Fig. 5.15 shows a continuous-fiber reinforced microstructure with the 
same properties as for Fig. 5.14. We vary the resolution from 64 x 64 
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to 1024 x 1024 pixels. Similar to scenario 1, see Fig. 5.14a, we apply a 
macroscopic strain £22 for all resolutions, so that both damage cases, i. e., 
dilatation and distortion, are being activated. The strain is successively 
increased from 0% to 5% within 50 equidistant loading steps. 

The resulting distribution of the predominant damage variable sph for 
damage due to dilatation is shown in Fig. 5.15. We observe that areas of 
low and high damage level are captured also for low resolution, but there 
are slight differences in the achieved damage level. Also, as expected, 
localization behavior is not evident. 


(a) 64? (b) 128? (c) 256° 


Asph 
0.15 
0.10 
0.05 
0.00 


(d) 512? (e) 1024? 


Figure 5.15: Model response for a continuous-fiber reinforced microstructure evaluated at 
five different resolutions 


To get a more qualitative insight, the macroscopic stress-strain curves 
are shown in Fig. 5.16a. For a resolution of 64 x 64 pixels, the computed 
stresses are overestimated. For higher resolutions with 128 x 128 
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to 1024 x 1024 pixels, the differences are small. Investigating the relative 
deviations (535° — 51924) /a48?4 of the computed effective stress 22 
relative to the stress at a resolution of 1024 x 1024, see Fig. 5.16b, we 
observe that, for resolutions with 256 x 256 pixels and higher, the 
deviations are below 1%. 
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Figure 5.16: Resolution study for the continuous-fiber reinforced microstructure 


The iteration counts and timings are collected in Tab. 5.2. The total outer 
iterations (including Newton and CG iterations) for all loading steps 
vary in a narrow window around approximately 2950 for all resolutions 
considered. To give a comparable standard for the inner (material) 
iterations, we compute the average number of inner iterations over 
all voxels and subsequently take the maximum value over all loading 
steps. A value of about 2.8 inner iterations per voxel, irrespective of 
the resolution, indicates quadratic convergence of the Newton method. 
Both, the overall timing for computing all inner iterations, as well as the 
total timing, increase roughly with the degrees of freedom. 
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Based on this resolution study, a resolution of 256 x 256 pixels represents 
a compromise between an accurate prediction and a short runtime, see 
Sec. 5.4.3. 


Table 5.2: Iterations and timings for the conducted resolution study 


resolution iterations u | timings i 
| total outer max. average inner | tinner tiotal 
64? 2953 2.774 2s 9s 
128? 2897 2.814 6s 26 s 
256? 3012 2.843 22 s 108 s 
5122 3023 2.853 88 s 524 s 
1024? 2 987 2.858 352s 1793s 


5.4.4 Tensile tests of bone specimens with different fiber 
bundle orientations 


As a first macroscopic example, we compute the predicted damage 
evolution for a tensile bone specimen (see Fig. 5.17) with aligned E- 
glass fibers in a surrounding UPPH matrix. We assume the fibers to be 
arranged as fiber bundles with a volume fraction of 50%, as discussed 
in Gorthofer et al. (2020), see also Sec. 4.3.2. The resulting (effective) 
transversely isotropic stiffness parameters are listed in Tab. 5.3. 


Table 5.3: Material parameters for the transversely isotropic SMC composite fiber bundles 


EL Er Gur VLT VIT 
37.73 GPa 10.33 GPa 3.64GPa 0.292 0.477 
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We model the tensile bone specimen in the commercial FE software 
ABAQUS/Standard (Dassault Systèmes, 2014) by a C++ user material 
subroutine, as stated in Sec. 5.4.1. The specimen is subjected to a 
displacement-controlled loading in e;-direction up to a strain of 5%, 
as shown in Fig. 5.17a. We vary the orientation of the fiber bundles 
from 0° to 90° in the eı-e>-plane (see Fig. 5.17b) and evaluate the 
stiffness evolution due to damage. We consider all four Puck cases 
with damage-hardening parameters as listed in Tab. 5.1. The hardening 
parameter for case I is set to 230 MPa to delay the respective damage 
evolution. We exploit the symmetry of the tensile bone specimen and 
model only one eighth of the part, as shown in Fig. 5.17a. The average 
runtime for 11375 fully integrated continuum eight-node brick elements 
(C3D8) and 100 steps is about 120 s to 180 s on 6 threads. 


Fibers 


+ is modelled 
8 


Loading in e,-direction 


(a) Loading and dimensions in mm (b) Considered orientations in discussed analysis 
(see Fig. 5.18 - Fig. 5.21, cases (a) - (d) each) 


Figure 5.17: Tensile bone specimen with different material orientations 


We compute the direction-dependent Young’s moduli via equation (5.76) 
and normalized each w.r.t. the initial Young’s modulus in the cor- 
responding direction, similar to Sec. 5.4.2. This provides a scalar 
value E/E between 1 and 0 for each evaluated direction, with E/E = 1 
if the material (point) is sound and E/E» = 0 for a totally damaged 
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material (point). We plot the resulting distributions of the normalized 
Young’s moduli for the global coordinate system directions (e1, e2 and 
e3) in Fig. 5.18 - Fig. 5.20, respectively. 

In Fig. 5.18 we see the reduction of the normalized Young’s modulus in 
loading direction e,. The damage in bundle direction is proportional 
to the amount of fiber bundles that are aligned to the loading direction. 
As shown in Fig. 5.7 (see Sec. 5.4.2) damage in bundle direction results 
from stress in the respective direction based on Puck case I. If all bundles 
are aligned perpendicular to the loading direction (see Fig. 5.18d) the 
material is damaged due to Puck cases II and III. Based on the considered 
material parameters, the reduction of the normalized Young’s modulus 
is less pronounced for these cases compared to case I. Fig. 5.18b and 
Fig. 5.18c, with bundles being rotated 30° and 60° about the e3-axis, 
exhibit a combined damage evolution due to cases I, II and III. The 
resulting normalized Young’s modulus in e;-direction decreases with 
a decreasing influence of case I and subsequently increases again with 
increasing influence of cases IT and MI. 


E/E 
1.00 

| 0.50 

0.00 


(a) 0° rotation (b) 30° rotation (c) 60° rotation (d) 90° rotation 


Figure 5.18: Normalized Young’s modulus in eı-direction for different fiber bundle 
orientations 


The normalized Young’s modulus in e2-direction is shown in Fig. 5.19. 
We see a normalized Young’s modulus of value 1 for fiber bundles which 
are aligned to the e,-direction (see Fig. 5.19a) and for a rotation of 90° 
(see Fig. 5.19d). Either Puck case I or cases II and III are active in loading 
direction e; (but no combination!) and, therefore, the normalized 
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Young’s modulus in e2-direction for these two bundle orientations 
remains unaffected. If the bundles are aligned in 30° or 60° about 
the eı-axis, a combination of cases I - III will be active, leading to 
a reduction of the normalized Young’s modulus in e2-direction (see 
Fig. 5.19b and Fig. 5.19c). 


E/Eo 


NS ONIN 


(a) 0° rotation (b) 30° rotation (c) 60° rotation (d) 90° rotation 


Figure 5.19: Normalized Young’s modulus in e2-direction for different fiber bundle 
orientations 


The normalized Young’s modulus evaluated in e3-direction is shown 
in Fig. 5.20. Similar to the e2-direction, we observe no damaging in 
the considered direction, provided the bundles are aligned in loading 
direction eı, as only Puck case I is active. The corresponding normalized 
stiffness has the value 1 in Fig. 5.20a. Due to the definition of the 
Puck cases, involving an averaging procedure, stiffness degradation 
perpendicular to the fiber directions (cases II and III) is interlinked, 
as discussed in Sec. 5.3.4 and Sec. 5.4.2. Consequently, damage in e3- 
direction is proportional to the rotation angle of the bundle orientation, 
see Fig. 5.20b - Fig. 5.20d. This is a direct consequence of the increasing 
influence of Puck cases II and II. As the bundle orientation varies 
based on a rotation around the e3-axis, this direction always remains 
perpendicular to the fiber direction. Hence, the normalized Young’s 
moduli evaluated in directions e3 and e% differ, especially for a bundle 
orientation of 90°, compare Fig. 5.20d and Fig. 5.19d. 

In addition to evaluating the damaged stiffness w. r. t. the global axes 
(e1, ea and e3), we would like to analyze the stiffness degradation 
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E/E 
1.00 
0.50 
0.00 


(a) 0° rotation (b) 30° rotation (c) 60° rotation (d) 90° rotation 


Figure 5.20: Normalized Young’s modulus in e3-direction for different fiber bundle 
orientations 


w.r.t. the local axes, namely in bundle direction (0°, 30°, 60° and 
90°) and perpendicular to the bundle direction. Fig. 5.21 shows the 
distribution of the normalized Young’s modulus evaluated in the bundle 
direction. Provided the bundles are aligned to the loading direction eı, 
the local bundle axis and the global e,-axis coincide, as do the resulting 
normalized Young’s moduli, see Fig. 5.21a and Fig. 5.18a. 


The further we rotate the bundle direction away from the loading 
direction, the higher is the remaining normalized Young’s modulus 
in bundle direction. As damage in bundle direction is solely caused 
by Puck case I (see Sec. 5.4.2 and Sec. 5.4.2), the respective damage 
variables do not evolve if all fibers are aligned perpendicular to the 
loading direction, see Fig. 5.21d. 


SN NIN IN 


(a) 0° rotation (b) 30° rotation (c) 60° rotation (d) 90° rotation 


E/Eo 
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Figure 5.21: Normalized Young’s modulus in fiber direction for different fiber bundle 
orientations 
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For the setup considered, the normalized Young’s modulus perpen- 
dicular to the bundle direction is the normalized Young’s modulus in 
e3-direction, as already discussed above. Hence, the distributions of the 
normalized stiffness perpendicular to the fiber direction are shown in 
Fig. 5.20. If all bundles are aligned to the loading direction (see Fig. 5.20a), 
no damage perpendicular to the bundles will evolve, in contrast to 
Fig. 5.21. The normalized Young’s modulus perpendicular to the bundle 
direction increases with increasing orientation angle. If all bundles 
are aligned perpendicular to the loading direction (see Fig. 5.20d), the 
corresponding damage perpendicular to the bundle direction is highest, 
which is a direct consequence of Puck cases IT and Il. 


5.4.5 A plain-weave composite under shear loading 


AG12/G12 
20% 
60% 
100 % 


(a) Voxelized weave microstructure with (b) Relative reduction (5.81) of shear mod- 
four tows ulus G12 


Figure 5.22: Microstructure and predicted relative reduction of the shear modulus G12 in 
a plaine weave composite 


Last but not least, we demonstrate the utility of our model framework 
for modeling anisotropic damage evolution in a woven fiber reinforced 
composite. Simon et al. (2017) investigated the mechanical behavior of 
a plain-weave composite manufactured from continuous carbon fibers 
reinforcing an epoxy matrix resin, see Fig. 5.22a. The carbon fibers are 


130 


5.4 Computational investigations 


aligned unidirectionally in fiber tows that are regularly interwoven. As 
each of these tows consists of thousands of carbon fibers, it is customary 
to work with a multiscale scheme that considers three different scales: 
the macroscopic scale is large compared to the woven unit cell, see 
Fig. 5.22a, which constitutes the mesoscale. Within the latter, the tows 
are considered homogeneous and anisotropic. On the microscale, the 
tows get resolved in terms of continuous carbon fibers in an epoxy resin. 


The linear elastic moduli of the considered materials are listed in Tab. 5.4. 


Table 5.4: Elastic moduli of matrix, fibers and tows (Simon et al., 2017, Tab. 1 & 2) 


Young’s modulus | shear modulus | Poisson’s ratio 
in GPa in GPa 
epoxy | E = 3 G = 109 |v = 0.38 
n EL = 290 Gir = 20 VIT = 0.2 
ae a Gr = 9 |m = 01 
EL = 144 Gir = 2.58 YLT = 0.29 
tows |m = 74 [Gr = 191 | yp = 0.39 


These comprise the isotropic epoxy matrix and the transversely isotropic 
carbon fibers. The transversely isotropic engineering constants for the 
tows were obtained by linear elastic homogenization. Please note that the 
subscripts "L" and "T" refer to longitudinal and transverse, respectively. 
Based on earlier work (Stier et al., 2015; Bednarcyk et al., 2015; 2014), 
Simon et al. (2017) presented a regularized orthotropic continuum 
damage-model based on the framework developed by Barbero and 
co-workers (Barbero and Lonetti, 2002; Lonetti et al., 2003; 2004)) and 
concisely summarized in his book (Barbero, 2013). More precisely, their 
strategy takes the orthotropic engineering constants as the point of 
departure, and models their degradation on an individual basis in terms 


131 


5 A convex anisotropic damage model based on the compliance tensor 


of associated scalar damage variables. Based on the associated driving 
forces, damage surfaces are defined, together with appropriate kinetic 
laws. 

Expressing the dependence of the stiffness tensor on the orthotropic 
engineering constants is most easily realized in terms of the compliance 
tensor, the approach of Simon et al. (2017) appears superficially similar 
to our approach. However, we do not fix the damage variables a priori. 
Rather, they emerge naturally in our framework based on the chosen 
extraction tensors and damage-activation functions. 

In this paragraph, we demonstrate that our model is capable of 
reproducing the damage behavior upon quasi-static loading of the weave 
composite. The protocol we present is straightforward and proceeds 
step by step. As a first step, we introduce a number of extraction tensors 
which capture elementary damage cases evoked by pure normal- and 
shear-loading scenarios. The tensors extract the associated normal 
and shear stress components from the applied stress state ø. These 
extraction tensors are uncoupled. Hence, damage in a certain direction 
is solely driven by the associated loading case, e.g., normal damage 
in eı-direction due to normal loading in eı-direction, 


ii = est, 2. ey, 33 = en (5.78) 
23 = (€2 Qs €3)°”, 13 = (€1 ®s e3)°”, 12 = (eı ®s e2)” 
(5.79) 


Combining suitable damage-activation functions based on these ex- 
traction tensors permit modeling a wide range of damage-evolution 
predictions. In particular, they enable us to describe the stiffness 
reduction for the scenario considered by Simon et al. (2017). 

We first capture the damage evolution in a neat epoxy sample under 
non-monotonic uniaxial loading and choose an extraction tensor of 
type (5.78);. Subsequently, we account for the damage onset due to shear 
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loading by using a second damage-activation function in combination 
with an extraction tensor of type (5.79)3. The identified parameters for 
the epoxy damage-model are summarized in Tab. 5.5, which were chosen 
to reproduce the results of Simon et al. (2017) best. 

Furthermore, we employ a number of damage-activation functions 
and suitable extraction tensors to capture the damage evolution in the 
fiber tows. We fix the the longitudinal tow direction to correspond to 
the local e;-direction. The response to shear loading in longitudinal 
and transverse directions is best described with extraction tensors of 
the forms (5.79); and (5.79)3. As the reduction of the two orthotropic 
Young’s moduli in the transverse plane (and hence the associated dam- 
age evolutions) is not identical, we introduce an additional extraction 
tensor 


= +62", (5.80) 


which is supplemented by a fourth damage-activation function with an 
extraction tensor that drives damage in e2-direction (5.78)2 only. Tab. 5.5 
comprises a complete list of the identified damage parameters. 


Table 5.5: Extraction tensors and identified damage parameters to capture the mechanical 
behavior of epoxy and tow 


extraction tensor oj9inMPa HinMPa m 
(5.78), 0.5 180 0.47 
epoxy (5.79); 101 41 0.97 
(5.79), 8 275 0.3 
ee (5.79)3 10 245 0.3 
(5.78) 79 30 0.9 
(5.80) 200 630 0.4 


The listed extraction tensors and damage parameters at hand allow us 
to reproduce the structural behavior of, both, the neat epoxy and the 
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tows, the latter in terms of stress-strain curves and the reduction of the 
orthotropic engineering constants. The corresponding results are shown 
in Fig. 5.23. 
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Figure 5.23: Comparison of predicted stress-strain curve and reductions of the orthotropic 


engineering constants based on introduced extraction tensors (see Tab. 5.5). Our model 
predictions are dashed in (b)-(d) 
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With the introduced extraction tensors and damage parameters at hand, 
we are able to reproduce the experimental results obtained for the 
neat epoxy resin, as well as the predictions computed by Simon et al. 
(2017) quite accurately, see Fig. 5.23a. The decrease in the individual 
orthotropic engineering-constants are shown in Fig. 5.23b to Fig. 5.23d, 
where dashed lines correspond to our model and solid lines refer to the 
references (Simon et al., 2017; Bednarcyk et al., 2015). 


For all considered loading cases, our proposed framework makes 
it simple to account for those engineering constants which remain 
unaffected during the loading. Fig. 5.23b shows that the reduction 
of Young’s modulus E33 is predicted correctly up to an applied strain 
of 10%. The shear moduli are predicted quite well up to a strain level 
of 5%, see Fig. 5.23c and Fig. 5.23d. Only the evolution of Young’s 
modulus E22 shows some deviations beyond a strain level of 3%, as 
our model is unable to capture the slope of reduction accurately enough 
for this case. To sum up, the proposed framework permits reproducing 
the reduction of all affected engineering constants accurately for small 
deformations of up to 3%. 

For the identified parameter set listed in Tab. 5.5, we investigate the 
plain-weave composite that was so 
also studied by Simon et al. (2017), 
see Fig. 5.22a. The mesostruc- 
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regular grid with 512 x 512 x 56 strain curves for plain weave under shear 
voxels. loading 
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Just as Simon et al. (2017), we investigate the longitudinal shear response 
of the plain-weave cell. For this purpose, we analyze the effective stress- 
strain curve, as well as the reduction of the shear modulus G12 measured 
in terms of the relative error 


AG = G9 Cig (5.81) 


w.r.t. the initial, undamaged shear modulus G{,. The predicted full- 
field distribution of this relative reduction is shown in Fig. 5.22b, and 
coincides well with the results presented by Simon et al. (2017), Fig. 9. 
Moreover, the effective stress-strain curves of the plain-weave composite 
subjected to longitudinal shear match well for the entire loading regime, 
see Fig. 5.24. 


5.5 Conclusions 


In this article, a generalized standard material (GSM) model for 
anisotropic damage evolution based on the compliance tensor as the 
primary damage variable was developed. Based on the insight that the 
Hookean elastic energy density, considered as a function of the elastic 
strain and the compliance tensor, is a convex function of both arguments, 
a convex framework for quasi-static damage evolution was established, 
preventing damage localization intrinsically. Indeed, by choosing the 
energy (density) related to the progressive degradation of the material 
appropriately, the condensed incremental potential (Miehe, 2002) is 
strictly convex and of superlinear growth, which prevents localization 
for such a model. Of course, working with a softening-type energy for 
the damage-surface variables is also possible, and should be studied 
more closely in subsequent work. 

Section 5.2 is organized to emphasize the modular fashion that the 
compliance-based damage model is built up. The model might be 
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extended in subsequent work, for instance accounting for strain-rate 
sensitivity within the model. For an overview of the assumptions 
leading to specific specializations of the evolution of internal variables, 
we refer to the overview in Appendix B. The modeling framework is 
general enough to incorporate coupling to other inelastic models, such as 
plasticity or viscoplasticity (Rousselier, 1979; McDowell, 2008), entirely 
within the proposed framework. Also, due to its inherent stability, an 
extension to fatigue damage, as observed for certain fiber reinforced 
polymers (Sauer and Richardson, 1980; Bartkowiak et al., 2019; 2020), 
appears promising (Magino et al., 2022). 


The modular character of the model was exemplified by specific 
extraction tensors motivated by Puck’s anisotropic failure criteria (Puck 
and Schtirmann, 2002; Knops, 2008), see Sec. 5.3. With these ingredients 
at hand, we demonstrated the model’s capabilities of developing 
complex anisotropic stiffness states, see Sec. 5.4, not restricted a priori by 
a specific degree of (an)isotropy of the stiffness tensor, emphasizing that 
the model is capable of handling any initial stiffness. We also showed the 
model’s capabilities on meso and volume-element scale, based upon a 
straightforward numerical treatment. With these achievements at hand, 
accounting for additional failure criteria (Kaddour et al., 2004; Fritzsche 
et al., 2008; Bouhala et al., 2013) or coupling the model to phase-field 
fracture models (Miehe et al., 2010; Schneider et al., 2016a; Gerasimov 
and De Lorenzis, 2019) appears possible. 

Returning to our original motivation, i.e., modeling anisotropic damage 
of SMC composite materials, requires incorporating the presented mod- 
eling framework into a three-scale homogenization scheme (Anagnostou 
et al., 2018). The underlying fiber bundle mesostructure (Schöttl et al., 
2021b; Meyer et al., 2020; Dumont et al., 2007) has to be accounted for, 
and the model parameters have to be fitted to experimental data. For 
the latter purpose, a convenient experimental program is necessary 
(Schemmann et al., 2018c). 
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From a mathematical perspective, a thorough mathematical analysis of 
our model is desirable, whereas an extension to tension-compression 
asymmetry appears imperative in order to model load reversals. 
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A computational multiscale 
model for anisotropic failure of 
SMC composites’ 


6.1 Introduction 


6.1.1 Research contributions 


We present a holistic multiscale approach for constructing anisotropic 
criteria describing the macroscopic failure of SMC composites. The 
approach is based on full-field anisotropic damage evolution on the 
microscale, and the failure criteria emerge naturally by accounting for 
the evolving microscopic damage. 

Based on findings of scientific studies, we evaluate experimental 
investigations regarding distributions of ultimate strength and stiffness 
reduction of sheet molding compound composite specimens. We 
summarize the characteristic microstructure, recapture micro-computed 
tomography analysis and the identification of linear elastic properties. 
We use our anisotropic damage model introduced in Chapter 5 on 
the microscale. The model directly operates on the compliance tensor, 
captures matrix and bundle damage via dedicated stress-based damage 


1 This chapter is based on the publication “A computational multiscale model for 
anisotropic failure of SMC composites” (Görthofer et al., 2022a) 
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criteria and allows for a comparison of simulation and experimental 
results. To identify the damage material-parameters used in the non- 
linear and time-consuming full-field simulations, we utilize Bayesian 
optimization with Gaussian processes. 

We validate the full-field predictions on the microscale and subsequently 
identify macroscopic failure criteria based on distributions taken from 
experimental findings. We propose failure surfaces in stress space and 
stiffness-reduction triggered failure surfaces to cover both a structural 
analysis and a design process perspectives. 


6.1.2 Chapter structure 


In this chapter, we introduce a methodology for identifying anisotropic 
failure surfaces to be used in component-scale simulations based on 
computational micromechanics. The presented multiscale approach 
allows for a scale-transition of full-field damage evolution within 
matrix and fiber bundles to an experimentally supported identification 
of application-based failure criteria for the SMC composite on the 
macroscale. 

With the desire to apply SMC composites as load-bearing structural 
components in mind, we start this chapter by analyzing experimental 
investigations as well as simulation-based studies published by different 
research groups, cf., Fitoussi et al. (1996; 1998); Anagnostou et al. (2018), 
Chen et al. (2018b); Li et al. (2018a) and the International Research 
Training Group “Integrated engineering of continuous-discontinuous 
long fiber reinforced polymer structures” (GRK 2078) (Kehrer et al., 
2018; Trauth, 2020; Schemmann et al., 2018a; Schöttl et al., 2021a). We 
give a brief overview of the considered SMC composite and a set of 
representative experimental investigations on neat UPPH and SMC 
composite specimens in Sec. 6.2. An analysis of wCT-based in-situ 
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experiments offers further insight into the behavior of SMC composites 
during loading. 

To represent the anisotropic elasto-damageable behavior of SMC com- 
posites (Fitoussi et al., 2005; Ogihara and Koyanagi, 2010; Trauth et al., 
2017a), we choose a modular anisotropic damage model which is 
formulated in the setting of generalized standard materials (Görthofer 
et al., 2022b) and utilizes a convex dissipation potential (Halphen 
and Nguyen, 1975; Hansen and Schreyer, 1994), see Chapter 5. This 
framework relies upon modular damage-activation functions and 
stress extraction tensors, fulfills Wulfinghoff’s damage growth criterion 
(Wulfinghoff et al., 2017) and ensures a well-posed model precluding 
localization. We provide an overview of the model and the associated 
set of parameters in Sec. 6.3. To capture the anisotropic damage in the 
SMC composite fiber bundles, we develop a set of extraction tensors 
motivated by Puck’s theory (Puck and Schiirmann, 2002; Knops, 2008) 
accounting for maximum stress states. To model the UPPH matrix 
behavior, we introduce an extraction tensor accounting for damage due 
to dilatation. 

We identify all associated damage parameters via a Bayesian optimiza- 
tion approach with Gaussian regression, as presented in Sec. 6.4. Using 
an anisotropic kernel function and parallel executions of FFT-based 
full-field computational homogenization (Schneider, 2021) simulations, 
we are able to unambiguously identify damage parameter sets for the 
SMC composite in Sec. 6.5. 

In Sec. 6.6, we analyze our full-field damage predictions and compare 
the results with uCT scans provided by Schöttl et al. (2020); Schöttl et al. 
(2021a). Furthermore, we compute three-dimensional, effective failure 
surfaces for the SMC composite, using failure criteria based on failure 
distributions evaluated in the experimental investigations, see Sec. 6.2. 
Hence, we provide both full-field damage evolution in SMC composites 
on the microscale as well as macroscopic failure surfaces for different 
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failure criteria taking the anisotropic and heterogeneous microstructure 


of SMC composites into account. These micromechanics-based failure 


surfaces are formulated in stress space to be used in a structural analysis 
(Gorthofer et al., 2019b), as well as via residual stiffnesses to be used 


in a design (optimization) process (Revfi et al., 2021). Our presented 


multiscale approach helps paving the way for an application-based 


anisotropic damage and failure evaluation of SMC composites on 


component scale level. 


6.2 Experimental investigations 


6.2.1 Neat UPPH 


Within the IRTG 2078 neat UPPH 
specimens could be manufactured 
by M. Bartkowiak (Institute for Ap- 
plied Materials — Materials Science 
and Engineering), based on prelim- 
inary research by Kehrer, Trauth 
and co-workers (Kehrer et al., 2018; 
Kehrer, 2019; Trauth, 2020). The 
formation of bubbles within the neat 
UPPH, which lead to local stress 
concentrations and hence foster pre- 
mature failure, could be minimized, 
allowing for a meaningful investi- 
gation of said material not only in 
the elastic, but now also in the non- 
linear regime. 


60 


a 
& 
Z aob 4 
b 
a or | 
0 1 L L L 
0 0.5 1 1.5 2 
Strain £ in % 
— 0° specimens — 90° specimens 


Figure 6.1: Neat UPPH bone specimens 
loaded in different directions (with cour- 
tesy of M. Bartkowiak and A. Trauth, 
Institute for Applied Materials — Materials 
Science and Engineering) in analogy to 
investigations by Kehrer et al. (2018) 


Specimens oriented in different directions were cut from manufactured 


neat UPPH sheets and tested uniaxially. The resulting stress-strain 
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curves, as well as the shape of the specimens, are shown in Fig. 6.1. 
The measured Young’s moduli, according to DIN EN ISO 527-4 
(Technisches Komitee CEN/TC 249, 1997), coincide with the Young’s 
moduli measured by Kehrer et al. (2018), as shown in the histogram 


Fig. 6.2a. 


The computed mean value for the considered experiments is 3.38 GPa 
with a standard deviation of 0.06 GPa. Furthermore, the non-linear 


regime up to failure is visible. As the structural behavior of the different 


specimens cannot be distinguished, it is reasonable to use an isotropic 


90° specimens 


model for UPPH. 
10° specimens 
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Figure 6.2: Analysis of the neat UPPH bone specimens shown in Fig. 6.1 
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The ultimate strength, i.e., the stress at failure, has an average value 
of 52.6 MPa with a standard deviation of 5.5 MPa, see Fig. 6.2b. The 
bandwidth of values for the ultimate strength, also accounting for the 
different orientations, are highlighted as horizontal coloring in Fig. 6.1. 
Assuming a pure elasto-damageable behavior, we determined the 
Young’s moduli at failure as a measure of stiffness reduction (or damage 
evolution) via Ao /Aeé, see Fig. 6.2c. The mean of all Young’s moduli at 
failure is 3.01 GPa with a standard deviation of 0.08 GPa. Relating the 
Young’s moduli at failure to the corresponding initial Young’s moduli we 
get the stiffness reduction shown in Fig. 6.2d. On average, this reduction 
is 10.67 % with a standard deviation of 3.19 %. 


6.2.2 SMC composite 


For our SMC composite, we follow a similar procedure to the inves- 
tigations on neat UPPH specimens. Based on previously conducted 
experimental investigations in terms of dynamic mechanical analysis 
(Kehrer et al., 2018; 2020), as well as loading and unloading scenarios, 
conceivably coupled to acoustic emission analysis (Trauth et al., 2017a; 
2021; Trauth, 2020), we know damage to be the predominant mechanism 
in the considered high performance SMC composite. Hence, we neglect 
viscous and plastic effects and focus on the analysis of uniaxial tensile 
tests to evaluate the onset of damage. Therefore, specimens with 
different orientations were cut from manufactured SMC composite 
sheets with a fiber volume fraction of 25 %. These sheets were produced 
in a two-dimensional compression molding flow-process where the 
bi-staged sheets are placed in the middle of the tool. Four different 
specimen geometries were used, see Fig. 6.3 - two rectangular specimens 
(types A and B) and two bone specimens (types C and D) with widths 
of 15 mm and 30 mm, respectively. 
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Figure 6.3: SMC composite specimens of different geometry types, loaded in different 
directions, investigated by Trauth, Bartkowiak and co-workers (Trauth et al., 2017a; Trauth, 
2020) 


Fig. 6.4a shows the measured Young’s moduli corresponding to the 
stress-strain curves in Fig. 6.3. These Young’s moduli have a mean 
of 11.85 GPa with a standard deviation of 1.92 GPa, and confirm the 
moduli measured by Trauth et al. (2017a) and Kehrer et al. (2018). The 
measured ultimate strength has a mean of 144.78 MPa and a standard 
deviation of 15.15 MPa, as shown in Fig. 6.4b. The total bandwidth 
of the ultimate strength is highlighted in Fig. 6.3. The histogram and 
normal distribution of the damaged (reduced) Young’s moduli at failure 
are shown in Fig. 6.4c. On average, the damaged Young’s moduli 
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take a value of 8.85 GPa with a standard deviation of 0.94 GPa. The 
average stiffness reduction at failure is 23.67 % with a standard deviation 
of 12.79 %, as shown in Fig. 6.4d. 
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Figure 6.4: Analysis of the SMC composite specimens shown in Fig. 6.3 


6.2.3 In-situ uCT scan analysis 


In addition to pre-processing uCT scans on unscathed specimens, also 
in-situ CT scans were performed by Schöttl et al. (2020); Schöttl 
et al. (2021a). Snapshots of a specimen loaded vertically in such an 


in-situ experiment, aiming mainly at damage evolution and microcrack 
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segmentation, are shown in Fig. 6.5. The characteristic fiber bundle 
microstructure is clearly visible in the pCT scan shown in Fig. 6.5a. 


(a) Initial (b) Loading step 1 (c) Loading step 3 (d) Loading step 4 


Figure 6.5: Snapshots of an in-situ SMC composite specimen for different loading steps in 
vertical direction showing the evolution of microcracks, as investigated by Schöttl et al. 
(2020) and Schöttl et al. (2021a). All shown snapshots are published under CC BY 4.0 
license in Schöttl et al. (2021a). 


Cyclic loading with increased loading steps and in-between holding 
times for CT scan analysis evoke a relatively stable damage evolution 
within the material due to microcrack propagation and hence allow for 
damage detection via scanning (Schöttl et al., 2020). The investigations 
show matrix damage perpendicular to the loading direction to be 
the primary damage mechanism (Schöttl et al., 2020; Trauth et al., 
2017a). Microcracks initiate at rims of matrix-rich areas — mainly 
at the edge or at fiber bundle interfaces — and progress through the 
material, diverted by bundles (Schöttl et al., 2021a). Depending on the 
orientation of individual bundles and hence the associated local stress 
state, microcracks either follow the principal bundle direction or run 
through it. Increasing the loading level leads to a higher microcrack 
density and a wider spread, until eventually total macroscopic failure 
is inevitable, see Fig. 6.5c and Fig. 6.5d. Fiber or bundle breakage is 
hardly observed and partially occurs at the point of macroscopic failure 
(Meraghni and Benzeggagh, 1995; Trauth et al., 2017b). 
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6.3 A modular framework to describe 
anisotropic damage based on extraction 
tensors 


6.3.1 A compliance-based anisotropic damage model 


To describe damage evolution within the SMC composite, we implement 
a corresponding model (Görthofer et al., 2022b) formulated in the 
framework of generalized standard materials (GSM) (Halphen and 
Nguyen, 1975). The modular nature of the model allows for an effortless 
employment to describe both damage within the matrix and the fiber 
bundles, based on associated extraction tensors. In the following, we 
briefly recall the most important aspects of the damage model (Gorthofer 
et al., 2022b) as presented in Chapter 5. 

The first potential in the GSM framework is the free energy in terms of a 
Hookean elastic energy and an energy related to the progressive damage 


w : Sym(d) x Sa x R” > R, 
M 
1 —1 H; mi+1 
(es, > 5e-S [e] +X , 4; , m;>Q. 


(6.1) 


The former is jointly convex in the strain e and the compli- 
ance S and infinitely often differentiable. The latter is intro- 
duced as a power-law hardening type in terms of M differ- 
ent damage variables q;, hardening parameters H; and expo- 
nents m; and hence is convex and continuously differentiable. The 
space of symmetric dx d infinitesimal strain tensors e is Sym(d) 
(with d being the dimension, i.e., d € [2,3]), the compliance tensors 
are Se S4= {S € Sym(Sym(d)) |T -S[r] > 0 for all r € Sym(d)\{O}}, 
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and the damage variables q € R™ are M scalars in real space. A more 
detailed discussion can be found in (Görthofer et al., 2022b). 


The second potential in the GSM framework is the force potential 


0, @(T,8)<0 forall i=1,...,M, 


i (6.2) 
+oo, otherwise, 


which we introduce via M convex and continuously differentiable 
damage-activation functions 


Qi : Sym(Sym(d)) x R > R, 


(6.3) 
(T, 6;) = 2T - B? — of; + Hißi, i=1,...,M. 


Each damage-activation function is formulated in terms of the conjugate 
driving forces T € Sym(Sym(d)) for the compliance and ß; € R for the 
damage variables, respectively. Furthermore, a fourth-order extraction 


tensor B; € L(Sym(d)) and a damage-activation threshold oo; are 


introduced per damage-activation function ġ;. As the driving forces are 
given by their potential-based relations T = 30 ®oand b; = —Hiq;", 
we simplify the damage-activation functions ¢; in terms of the stress 7 


and the damage variables q; as 


fi: Sym(d) x R > R, 
(0,4) > |Bile] |? — 90; - Hi, i=1,...,M. 


ı 


(6.4) 


Biot’s dual equation (Borwein and Lewis, 2006) in the context of the 
introduced potentials (6.1) and (6.2) yields the evolution equations for 
the internal variables S and q, where the consistency parameters need 
to obey the classical Karush-Kuhn-Tucker (KKT) conditions (Karush, 
1939; Kuhn and Tucker, 1951). With a little work, we eliminate the 
consistency parameters and reformulate the equations, s. t. we compute 
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the compliance at time t as 


M (t) 
S(t) =S) +25- m 2 (6.5) 
i=1 


where Sp = S(0) is the initial compliance. The associated KKT conditions 
referring to the simplified damage-activation functions f; and damage 
variables q; read 


filo, qi) < 0, di 2 0, di file, qi) = 0, i= 1, we „M. (6.6) 


6.3.2 Influence of the model parameters 


The model is able to describe progressive, fully anisotropic damage of 
any (anisotropic) initial stiffness Cp = Sp‘ via any number of damage- 
activation functions (6.4). Each damage-activation function f; involves 
three parameters, a damage-activation threshold oo,;, a hardening 


parameter H; and a power-law exponent m;. As the squared norm of 
I? 


the extracted stress ||B; [a] ||? reaches of ;, damage evolution is triggered, 
and the damage variable q; increases. The damage evolution and 
hence the non-linear hardening-type stress-strain relation are governed 
by H; and m;. The influence of the parameters on the predicted 
stress-strain relations is shown in Fig. 6.6 as an example. Here, a single 
damage-activation function f to describe damage in loading direction, in 
analogy to Govindjee et al. (1995), is used. We choose the set of damage 
parameters as oo = 30 MPa, H = 130 MPa and m = 1 for reference. For 
each study shown in Fig. 6.6, we vary one of these parameters and keep 


the others constant. 

Increasing the damage-activation threshold ou retards the damage 
initiation and increases the elastic region, see Fig. 6.6a. The hardening 
parameter H and the power-law exponent m have an opposing influence 
on the damage region, as we observe when comparing Fig. 6.6b and 
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Fig. 6.6c. Whereas an increase of H increases the hardening-type 
behavior, an increase of m decreases the hardening. 
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Figure 6.6: Influence of the model parameters oo, H and m on the predicted stress-strain 
relation 


6.3.3 Extraction tensors accounting for dilatation and dis- 
tortion 


Spherical stress state 


To model damage induced by dilatation, i.e., due to spherical stress 
states, the corresponding extraction tensor B is the spherical projector P4, 
well-known from linear elasticity 


1 
=Pı=,I81. (6.7) 


Due to its characteristic P; [o] = o°, the extracted stresses of the defined 
damage-activation functions (6.4) read 


la] ||? = llo° ||. (6.8) 
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Deviatoric stress state 


In analogy to spherical stress states, we model damage induced by 


distortion, i. e., by deviatoric stress states, via an extraction tensor B of 
the form 


= Ps = I° — P. (6.9) 


Due to the characteristic of the deviatoric projector P2 [o] = o’, the 
corresponding extracted stresses are 


IB [ø] |? = lle’l?. (6.10) 


6.3.4 Puck-type extraction tensors accounting for maxi- 
mum stresses 


Basic idea 


Motivated by the basic stress states present in laminates, as investigated 
by Puck and co-workers (Puck and Schürmann, 2002; Knops, 2008), we 
define extraction tensors to account for damage due to normal as well 
as shear loading, both in fiber bundle direction and perpendicular to it. 
Let {e1, e2, e3} be a Cartesian coordinate system, where eı corresponds 
to the alignment direction of the fibers, see Fig. 6.7. Then the stress 
tensor o admits the block decomposition 


O11 | 012 013 


0O = 012 022 023 . (6.11) 


013 | 023 033 


The stress in fiber direction is given by o1ı, the lower right block 
describes the transverse stresses in the plane orthogonal to the fiber 
direction, and (012,013) are the remaining shear stresses in longitu- 
dinal fiber bundle direction. Given the basic stress state in a fiber 
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bundle (see Eq. (6.11)), the traction vec- Matrix 
tor tin an arbitrary direction n is t = on. 


Fibers 
To account for the maximum stresses in 


terms of a Puck-type setting, we need 
to maximize specific (extracted) stresses 
over certain directions using the pencil 
glide approach (Krawietz, 1981; 2001). a en o 
We distinguish four different loading tem {e1, e2, e3} 

scenarios that yield four associated ex- 


traction tensors. 


Normal loading in fiber direction 


The traction vector pointing in fiber direction eı is 
tı = (ey 3 oe) ei, (6.12) 


with the scalar stress 
ti = €] : O€ = 011. (6.13) 


This stress is unambiguously defined by the given bundle direction e1. 
The extraction tensor capturing normal loading in fiber direction is 


zen (6.14) 


Accordingly, the extracted stress reads 


IB [ø] |? = on. (6.15) 


In the context of Puck’s theory on laminates we apply this extraction 
tensor to capture damage evolution in fiber bundle direction. Nontheless, 
it is not limited to fiber bundles and can be generally applied to describe 
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damage due to normal loading, similarly to the introduced tensor by 
Govindjee et al. (1995). 


Normal loading perpendicular to fiber direction 


An arbitrary traction vector perpendicular to the fiber direction is given 
via the direction $? 5 k L eı, where S? = {ze R? | ||x|| = 1}, as 


to = (k-ok)k. (6.16) 
We need to maximize the normal stress pointing into direction k 
to = k- øk, (6.17) 
over all admissible directions 


toes = max k- ok. (6.18) 
ke, 


Hence, t#?** is the maximum principal stress of the stress state on the e1- 
plane. This stress state is described by 


Og = Oijjei ® ej, 1,9 = 2, 3, (6.19) 


with the corresponding principal stresses 


2 
og = i 33 | Ne - z) ge, (6.20) 


Following the classical convention for principal stresses Aa > Az, we get 


max — Yo, (6.21) 
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We derive the associated extraction tensor and the damage-activation 
function via the condition 


[or] |? = (em? (6.22) 


While this condition (6.22) does not allow for an explicit closed-form 
description of the associated extraction tensor, the extracted stresses are 


given via 


2 
ø]? = tn, (ee). (623) 


Based on the procedures and results for the alternate loading scenarios 
(see Sec. 6.3.4), as well as extraction tensors formulated in an average 
stress manner Görthofer et al. (2022b), we postulate a damage extraction 
tensor for normal loading perpendicular to the fiber direction of the form 


v2 2 V2 2 
= * (ef? + ef?) + (ef? — ef?) + (e2 8s e3)”. (6.24) 


The associated extracted stress therefore reads 


IB [ø] ||? = (5032 + 5033 + 6099033 + 4033) : (6.25) 


el = 


Shear loading perpendicular to fiber direction 


In order to evaluate the stress state due to shear loading perpendicular 
to the fiber direction, we need to consider the shear part of the traction 
vector, namely 

t-=ok-(k:ok)k. (6.26) 
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We decompose the shear traction vector into shear stresses perpendicular 
and parallel to the fiber direction 


t; = tri +t. (6.27) 
The shear traction vector parallel to the fiber direction is 
t- = (ei $ t,) eı = (eı . ok) €l; (6.28) 


and, analogously, the shear traction vector perpendicular to the fiber 
direction is 


t,, =ok-(k:ok)k-(e,: ok) eı 
Sass ~~ 


t, tai 
= (I-k®&k-e,®eı)ok 
= (m: ok) m, (6.29) 


where S? >m L k L e and m L e1. The sought scalar shear stress is 
given as 
t-ı =m-ok, (6.30) 


which we need to maximize over all admissible directions m L k L eı 


"= max m-ok. (6.31) 


Note that the directions m and k are not independent. For a given 
bundle direction eı, we express m as a function of k via the condition 
I=m®m+k®k+eı®eı. Using this relation we get the expression 


max — max y ||øk||? — (k - øk)? - (e,-ok)?, (6.32) 
ke, 
involving only a single unknown unit vector k. 
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Analogously to normal loading perpendicular to the fiber direction, 
t is the maximum shear stress on the plane defined by the bundle 
direction e1. In accordance with principal stresses computed for the case 
above, the maximum shear stress is 


2 
max L [A2 5 A3| _ € 5 z) pat. (6.33) 


The associated extraction tensor and the damage-activation function are 


computed via the condition 


! max 

IB [o] |? = ary’, (6.34) 

which is ensured, provided 
2 jo] + 2 7 783 (e3? — e3?) + 093 (e2 Qs es) (6.35) 
holds. Evaluating this condition yields an explicit form of the extraction 

tensor 
2 

= y? (= en" + (€2 8s es)”, (6.36) 


and the corresponding extracted stress 


IB [ø] |? = 


(039 + 033 — 2022033 + 4053) - (6.37) 


Ale 


Shear loading in fiber direction 


We account for the maximum stress due to shear loading in fiber 
direction via the corresponding shear traction vector (6.28) 


t-j = (eı : ak) eı = tre. (6.38) 
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We maximize the associated shear stress over all admissible direc- 
tions k l eı 


mex = max e- ok, (6.39) 
Le, 


i.e. we search the normal direction k that points into direction of the 


maximum stress (Krawietz, 1981; 2001). As the stress is symmetric 


o' = ø, the relation 


t =e1:ok=k-oe, (6.40) 


holds. We define a corresponding projector P; with the characteris- 
tics P,k = k and P} = P,as 


P,.=I-eı®eı. (6.41) 
Application of this projector to the stress state yields 
zal = k “oe, = Pk -Oei = k i Poe = 012€2 + 0713€3. (6.42) 


The shear stress is maximized whenever 


ı Poea 


k= SAn 
|Preeil 


(6.43) 


holds, yielding 


Poe 
max — —~____. Pk = ||P = 2 Tas 6.44 
|| [Picea || koe = |Proeill = y 712 + 013 (6.44) 


The essential relation 


[o] |? = Coe (6.45) 


holds, if 
! 


* [ø] = 012 (€1 Qs e2) + 013 (€1 Qs e3) - (6.46) 
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Evaluating this relation gives an explicit form for the associated 
extraction tensor 


B = (e1 8s €2)°” + (e1 8s es)”, (6.47) 
and the corresponding extracted stress 


IB [ø] ||? = of) + 013- (6.48) 


6.4 Bayesian optimization process 


The presented framework takes damage parameters into account (see 
Sec. 6.3.1 and Sec. 6.3.2), which need to be identified to describe the 
elasto-damageable behavior of our SMC composite properly. The 
model and its parameters serve as input for the FFT-based full-field 
homogenization computed with our in-house code homKIT (see Fig. 6.8 
and Sec. 6.5.1 for further information). The stress-strain response 
computed by homKIT is compared to the experimental results discussed 
in Sec. 6.2 via a corresponding error measure, yielding an optimization 
problem for the identification of proper damage parameters. 


homKIT 
= 
| sampling | p microstructure 


error measure 


experiments 


Figure 6.8: Simplified schematic workflow of parameter identification routine 
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We choose a Bayesian optimization approach (Mockus, 1989; 1994), 
as our problem at hand involves an objective function that is expen- 
sive to evaluate (each call involves a full-field homogenization on a 
three-dimensional voxel image) and derivatives that are not easily. 
Similar conditions hold for, e.g., ductile phase-field fracture (Noii 
et al., 2021). Bayesian optimization uses a surrogate model of the 
actual problem accounting for probabilities and uncertainties based 
on Bayesian statistics (Bolstad and Curran, 2016). 


We describe the problem in terms of an abstract cost function 
c:R>R, p> clp), (6.49) 


combining the full-field homogenization based on a parameter 
set p E€ R C RP, and the error measure. The feasible domain R for the 
parameter set p represents the parameter ranges and is frequently chosen 
as hyper-rectangle with p € R C R? : l; < pi < ui, i = 1,. .. ,b, for scalar 
lower and upper bounds l; and u;. We seek a minimizer p of the cost 
function 


in. 6.50 
c(p) on, (6.50) 


In general, we do not know the (continuous) cost function, 
but only a certain amount of discrete values for parameter 
sets pi. =(pi,---;Pr) ER”, which we assemble in a vector 
Cir = [c(p1),-..,c(pr)|' € R”. The cost function values for the 
remainder of parameter sets are uncertain, yielding an infinite number 
of possible cost functions. Hence, we use a surrogate cost function 
model based on Gaussian process regression (Williams and Rasmussen, 
2006) to account for the uncertainty within the Bayesian optimization 
procedure. We assume all possible cost function values to be distributed 
normally 

cp) ~ N(m(p), K(p,p')) forall pe R, (6.51) 
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with a mean m: R — Rand a kernel (covariance) k : R x R > R. The 
kernel k defines the correlations of two parameter sets p and p’. 

Condition (6.51) has to hold for all parameter sets p € R yielding a 
normal distribution of all cost functions. Therefore, given the already 
known parameter sets p., with associated cost function values cı:r, 
all possible (surrogate) cost functions have to be drawn according 
to the multivariate normal distribution (6.51) with the mean vector 


mi: = [m(p1),...,m(p,)]' € R” and the covariance matrix 
k(pı,pı) kes k(p1, Pr) 
Kinar = : = : € Mr, (6.52) 
k(pr, pı) sree k( pr, Pr) 


The multivariate normal distribution (6.51) has to hold for any new 
parameter set p,+1, as well. Hence, the joint distribution based on all 
parameter sets - known and new - (p1:,,Pr+1) follows a multivariate 
normal distribution that reads 


r K wlr K 7,7 
[GE ~N (0, 1 1 L:r,r+1 (6.53) 
c(Pr+1) ep act k(Pr+1, Pr+1) 
with the kernel evaluations K1:,.r+1 = [k(p1, Pr+1); - - - ; k(Pr,Pr+1)] and 


k(pr+1,Pr+1). For convenience we assume the mean to be zero mı.- = 0 
Frazier (2018). 


Consequently, the (conditional) posterior probability distribution of the 
cost function c(p,+1) for the new parameter set p,+1, given the already 
known c1:r, is computed as 


c(Pr+1) | Cir =N (™(pr41), Y (pr+1)) (6.54) 


with a posterior mean m(p,+1) and a posterior variance 07(p,41). Here, 
the Sherman-Morrison-Woodbury formula (Hager, 1989) can be utilized 
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(Williams and Rasmussen, 2006) and yields 


m(Pr+1) = Ken ee (6.55) 
0° (p41) = k(Pr+1, Pr+1) — Kleri Kirin erp (6.56) 


The posterior mean m(p,+1) is estimated on the given, known data c1:r 
weighted by the kernel. The posterior variance V (pr41) is given 
by the prior covariance which is corrected by a term that takes the 
known correlations into account. The posterior mean approximates the 
actual cost function and is updated and improved with each evaluated 
parameter set p,+ı. For an infinity number of evaluations r + oo, we 
recapture the actual (deterministic) cost function. 

As we know the different parameters of the different damage-activation 
functions to have different levels of relevance on the anisotropic behavior 
of our SMC composite, we use the Matern 5/2 kernel (Matérn, 2013) 


kipp) =% (1 +55) + 58 (wp) exp (v55 (p, p')) (6.57) 


with an anisotropic distance function ô (p, p') between parameter sets p 
and p’ 


b (ie n2 
dipp) = a Di (6.58) 


i=l 


The Matern 5/2 kernel (6.57) is not restricted to smooth cost functions. 
Each of the hyperparameters (a1,..., a) is determined to represent the 
relevance of its associated damage parameter. The hyperparameter 72 
controls the width of the kernel, i.e., the general level of correlation 
between two parameter sets. For a more detailed discussion on kernel 
functions, the reader is referred to the book of Williams and Rasmussen 
(2006). 
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For each evaluation, we want to choose a new parameter set p,+1 
that is closer to the sought minimum of the cost function. Hence, a 
promising parameter set p,+1 has a large improvement (c(p,.) — ¢(pr+1)) 
w.r.t. the currently best parameter set p, = argmin,—,.c(p;). As we do 
not know c(p;+1) beforehand, we take the expected value € of the 
improvement given all already known evaluations c1.,, and maximize it, 
yielding 

Pr+1 = argmax EI(p), (6.59) 


with the "expected improvement" (EI) acquisition function (Mockus et al., 
1978; Brochu et al., 2010) 


EI(p) =€ ((c(ps) z c(p)) | Cir) (6.60) 


and Macaulay brackets (-) = max (0, -). We evaluate (6.60) in closed form 
(Brochu et al., 2010) and arrive at 


Ac(p) Ac(p) 
o(p) o(p) 


with the probability and cumulative density functions PDF and CDF. 


EI) = (eto) +3) PDF (SER) — jac onr (SR) (661 


The expected difference is defined as Ac(p) = c(p..) — m(p) — &, includ- 
ing a shift-parameter € that helps controlling the exploration and 
exploitation trade-off (Brochu et al., 2010). This trade-off between 
choosing (new) parameter sets with a high uncertainty (high v(p)) 
or a high expected quality (high Ac(p)) is the elementary part of all 
acquisition functions. 

To initialize a certain number of cost function values for the optimization 
process, we use a Latin-Hypercube sampling approach (McKay et al., 
1979) for the parameter sets within our hyper-rectangular domain R. 
The Bayesian optimization process is stopped after a given amount of 
iterations and the parameter set p, yielding the smallest cost function 
value c(p..) is considered as the optimal choice. 
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6.5 Identification of damage parameters in 
SMC composites 


6.5.1 Computational setup 


We integrated the proposed damage model as a user-defined subroutine 
in our in-house OpenMP-parallel FFT-based computational homoge- 
nization code written in Python 3.7 with Cython extensions (Behnel 
et al., 2011) and FFTW (Frigo and Johnson, 2005) bindings, as described 
by Schneider (2018) (homKIT, see Fig. 6.8). We use the discretization 
by Moulinec and Suquet (1998) and a Newton-CG scheme to solve the 
ensuing nonlinear systems of equations (Gélébart and Mondon-Cancel, 
2013; Wicht et al., 2020). We implemented the Bayesian optimization 
procedure in combination with the kernel evaluations and the error 
measure in terms of a Python 3.7 code, incorporating GPyOpt (Paleyes 
et al., 2016) (fitKIT, see Fig. 6.8). We executed the optimization 
iterations in parallel, using a batch sampling as introduced by Thompson 
(1933). 

For the computations we used either 6 - 12 threads on a desktop 
computer with 32 GB RAM and an Intel i7-8700K CPU with 6 cores 
and a clock rate of 3.7 GHz, or a workstation with two AMD EPYC 7642 
with 48 physical cores each, enabled SMT and 1024 GB of DRAM. 


6.5.2 Definition of the parameter set to be identified 


The linear elastic properties of the UPPH matrix and fibers are listed 
in Tab. 6.1. To derive the properties of a fiber bundle, we implemented 
two different approaches, numerical full-field homogenization of a 
representative bundle (Moulinec and Suquet, 1998), as well as mean-field 
Mori-Tanaka homogenization (Mori and Tanaka, 1973; Duschlbauer et al., 
2003), see also Sec. 4.3.2. Both approaches lead to approximately the 
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same transversely isotropic stiffness properties that are listed in Tab. 6.1, 
where "L" and "T" refer to longitudinal and transverse, respectively. 
Measurements and comparisons of the effective elastic properties of 
SMC composites were performed by Trauth (2020). 


Table 6.1: Elastic material parameters for the considered SMC composite 


E-glass fibers UPPH matrix Fiber bundles 

(Kehrer et al., 2018) (Kehrer et al., 2018) 

Eriso = 72 GPa EM iso = 3.4 GPa Egy = 37.73 GPa 
Epr = 10.33 GPa 

VE iso = 0.22 M,iso = 0.385 VB,TT = 0.477 
VB,LT = 0.292 


GEiso = 29.51 GPa GM iso = 1.23 GPa GET = 3.58 GPa 
GBLr = 3.64 GPa 


To capture matrix and bundle damage, the dominant mechanisms in 
SMC composites (Fitoussi et al., 2005; Ogihara and Koyanagi, 2010; 
Trauth et al., 2017a), we we utilize the following three damage-activation 
functions (6.4). Matrix damage evolution due to normal loading in 
terms of dilatation is captured via a damage-activation function in 
combination with an extraction tensor of type (6.7). Bundle damage in 
transverse and longitudinal direction, accounting for normal and shear 
stresses, is described via damage-activation functions in combination 
with extraction tensors of types (6.24) and (6.47), respectively. 

Each damage-activation function comes with three parameters (oo, H, 
m), as discussed in Sec. 6.3.2. We know the power-law exponent m to 
have a contrary influence on the non-linear regime compared to the 
hardening parameter H. To describe the behavior of SMC composites, 
we assume a linear damage hardening evolution and hence set the 
power-law exponent to m = 1 for each damage-activation function. This 
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leaves us with six parameters to be identified adequately in order to fit 
the predicted behavior of our model to the experimental observations. 
A summary of the applied extraction tensors, the associated parameters 
and their allowed ranges for the optimization process is given in Tab. 6.2. 


Table 6.2: Extraction tensors and damage parameters used for modeling the behavior of 
the SMC composite 


| Type | coin MPa Hin MPa 

UPPH | Sec.6.33 (67) | [5,50] | 50,400] 
Sec. 6.3.4 (6.24) | [5,40] [100,700] 

Bundle | 5.06.34 (647) | [5.40] [50,400] 


6.5.3 Preliminary study on neat UPPH 


In a preliminary study, we fit the material model with solely a single 
damage-activation function to the experimental results on neat UPPH 
specimens as presented in Sec. 6.2.1. We use a representative experiment 
showing the mean behavior of all conducted experiments on UPPH, see 
Fig. 6.1. 

As known from previous investigations (Kehrer et al., 2018; Kehrer, 2019; 
Trauth, 2020), neat matrix specimens behave slightly different to the 
matrix within the composite. Hence, we apply an extraction tensor of 
type (6.14). Furthermore, the preliminary study does not replace the 
UPPH parameter identification in a coupled approach of matrix and 
bundles within the SMC composite. Nonetheless, a preliminary study 
helps to understand the general material behavior of UPPH and can also 
be considered a stand-alone damage parameter identification for the 
neat UPPH. 
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Utilizing the Bayesian optimization process discussed in Sec. 6.4, we run 
an optimization with 100 initial hypercube samplings and subsequent 
3000 evaluations, using 150 batches of 20 parallel executions each. We 
update the surrogate model and hence the acquisition function after 
the computation of each batch. The computations are performed on a 
single voxel with isotropic material parameters for UPPH as presented 
in Tab. 6.1. Based on knowledge gathered in trail computations, we 
set the parameter ranges for the optimization to oo € [0,30] MPa and 
H e [500,600] MPa. 


We use a least square ansatz for the cost function value regarding a 
parameter set p 


elp) = | < 5 (oz? — gsim)* (6.62) 


i=l 
comparing the difference of the experimentally measured and the 
predicted stress-strain curves. We evaluate the measured stress où? 
and the computed stress oŝ™ for N given strain values. 

The cost function values c(p) of a representative optimization run 
with its total 3100 iterations are shown in Fig. 6.9. The associated 
parameter sets p = {o0, H} are plotted in Fig. 6.10. The lighter lines 
show the computed results of all 3100 iterations, whereas the darker 
lines represent the best results per batch, i.e., the parameter set that 
leads to the minimum cost function value per batch. 

Within the first 100 iterations, the hypercube sampling is evident as 
possible parameter sets of the hyper-rectangle R are scanned in a 
regular fashion (see Fig. 6.10) and hence the cost function values vary 
strongly (see Fig. 6.9). For the next batch of 20 more iterations, the 
cost function varies around a value of approximately 0.06 MPa with an 
95 % confidence interval of about +0.01 MPa. After the next update of 
the optimization surrogate model, the cost function values converge to 


approximately 0.058 MPa for the remainder of the optimization run with 
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Figure 6.9: Values of the cost function c(p) over iterations 


an associated parameter set as shown in Fig. 6.10. Solely after about 1100 
iterations an alternative parameter set (way off the other sets promising 
to be an optimal choice) was tested and then discarded, yielding a 
short increase of the cost function. The choice of said parameter set 
is evoked by the batch sampling discussed below. In all following 
iterations, the 95 % confidence interval of the computed cost function 
value is of the order of 10~° MPa, i.e., we are very confident that our 
result is correct. The associated values of the parameter set converge to 
values of about oo = 12.07 MPa and H = 555.95 MPa, see Fig. 6.10a and 
Fig. 6.10b. 

The remaining variation of the parameter values and the associated 
gradual behavior of the cost function considering all iterations (see 
lighter lines in Fig. 6.9 and Fig. 6.10) is caused by the batch sampling 
(Thompson, 1933). Here, for each batch, the most likely best parameter 
set is chosen together with 19 other parameter sets sampled via a beta 
distribution, to ensure an appropriate balance of the exploration and 
exploitation trade-off (see Sec. 6.4). Hence, the immediate quality is to 
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Figure 6.10: Values of damage-activation threshold og and hardening parameter H over 
iterations 


be improved together with the aim to gather new information possibly 
improving a currently high uncertainty. 

If we only consider the best result per batch (see darker lines in 
Fig. 6.9 and Fig. 6.10), we find the mentioned variations to reduce 
significantly after the initial hypercube sampling and a few following 
iterations. Hence, we basically found the optimal parameter set with the 
corresponding minimal cost function value. 

To check convergence of our algorithm, we repeated the discussed 
optimization process in 100 runs with different initial samplings. All 
computed optimal cost function values and associated parameters 
coincide up to computational accuracy. Hence, we can reliably repeat 
the computation of the best parameter set and assume the values 
of oo = 12.07 MPa and H = 555.95 MPa (with a cost function value 
of 0.058 MPa, which is 0.11 % of the mean ultimate strength) to be correct. 
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A comparison of the experimental 
stress-strain curves as presented in 
Sec. 6.2.1, cf. Fig. 6.1, and the com- 80 
puted stress-strain curve is shown 
in Fig. 6.11. The applied damage- 
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framework, in combination with the 


an accurate prediction of the macro- 
scopic stress-strain relation of neat Figure 6.11: Comparison of effective stress- 
UPPH. strain curves 


6.5.4 Study on SMC composite microstructures 


Similarly to the pre-study conducted on neat UPPH in the previous 
section 6.5.3, we optimize the parameter set needed to describe proper 
damage evolution in the SMC composite (see Sec. 6.5.2). The allowed 
ranges for this parameter set during the optimization process are shown 
in Tab. 6.2. We chose these ranges based on trail studies to provide an 
optimization hyper-rectangular space that is as large as necessary but as 
small as possible. 

We ran our optimization with 300 initial hypercube samplings and 
subsequent 3000 iterations, computed in 250 batches with 12 parallel 
executions on 8 threads each. We update the surrogate model and hence 
the acquisition function after the computation of each batch. 

The underlying microstructure is a generated unit cell (Görthofer et al., 
2020), as shown in Fig. 6.12a, with a fiber volume fraction of 25% and a 
planar isotropic orientation, matching the conditions of the specimens 
discussed in Sec. 6.2.2. 
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ey 


(a) diag(0.5,0.5,0) - Planar (b) diag(0.7,0.3,0)- Higher (c) diag(1,0,0) - Unidirec- 
isotropic orientation alignment in e,-direction tional orientation in 
e.,-direction 


Figure 6.12: Example of generated unit cells. Coloring indicates bundle orientation. Matrix 
is hidden. 


Using the introduced microstructure generation approach (Görthofer 
et al., 2020), we transfer our results to SMC composite unit cells with 
different orientations, see Fig. 6.12 and Sec. 6.6. Therefore, we generated 
additional unit cells with fiber bundles oriented in e,-direction to a 
higher extent and unit cells with a pure unidirectional bundle alignment 
in e,-direction, see Fig. 6.12b and Fig. 6.12c. The associated fiber 
orientation tensors of second-order (Advani and Tucker, 1987) are given 
in the captions. 

The cost function values c(p) of the 300 initial samples and the best result 
per batch for the subsequent iterations are shown in Fig. 6.13. The lighter 
areas mark the 95 % confidence interval for each result. We observe a 
high variation in the cost function values for the 300 initial sampling 
iterations, as we raster the complete hyper-rectangular parameter space. 
In a sense, we maximize the exploitation and minimize the expectation 
in our discussed trade-off (see Sec. 6.4), to gain as much knowledge 
about the cost function behavior as possible during initialization. 

In the following 1400 iterations the variation of the cost function 
decreases, but nonetheless remains relatively high, see Fig. 6.13a and 
Fig. 6.13b. As the behavior of the cost function is governed by six 
parameters, the individual influences need to be explored, which is 
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Figure 6.13: Best value per batch of the cost function c(p) over iterations 


responsible for this prolonged variation with its quite large confidence 
interval. After about a total of 1700 iterations, a possible optimum is 
found with a cost function value of approximately 0.45 MPa and a small 
confidence interval of the order of +0.01 MPa, see Fig. 6.13b. For the 
remainder of the optimization process, the variation of the best cost 


function value per batch is negligibly small with an associated small 
confidence interval, see Fig. 6.13c. 


Table 6.3: Extraction tensors and identified SMC composite damage parameters 


| | ooin MPa sin% | HinMPa sin% 
UPPH | (6.7) 43.99 24 177.84 3 
(6.24) | 29.36 11 381.64 4 
Bundle | (car, | 2860 51 123.81 7 


Hence, we can assume the best overall cost function value of 
about 0.44 MPa (which is 0.33% of the mean ultimate strength) to be 
reasonably close to the actual optimum. The corresponding values of the 
best parameter set to describe damage evolution in our SMC composite 
are listed in Tab. 6.3. Furthermore, the sensitivity s of each parameter in 
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terms of its associated hyperparameter related to all hyperparameters is 
added. In general, the identified optimum of the cost function shows a 
higher sensitivity w. r. t. the damage-activation threshold oo compared 
to the hardening parameter H. 

A comparison of the experimentally 
measured stress-strain curve (see 
Sec. 6.2.2, cf. Fig. 6.3) and the com- 
puted stress-strain curve is given in 


Fig. 6.14. These curves, representing ey 


the macroscopic structural behav- 
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6.6 Identification of macroscopic failure sur- 
faces for SMC composites 


6.6.1 Full-field damage evolution 


The introduced damage cases and identified parameters do not only 
allow for the prediction of the macroscopic stress-strain behavior, but 
also help us to evaluate the evolution and full-field distribution of 
damage on the microscale. With damage in terms of microcrack 
evolution (Schöttl et al., 2020; Schottl et al., 2021a) in mind, we analyze 
the predicted damage in the UPPH matrix. The evolution of damage over 
different loading steps is shown in Fig. 6.15. In analogy to the discussed 
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experiments, we consider a slice of a planar isotropic microstructure that 
is loaded in vertical direction. Damage in the matrix mainly initiates 
at the rims of matrix rich regions at bundle edges and corners, see 
Fig. 6.15a and Fig. 6.15b. With increasing loading level, on the one 
hand the damage level of existing damaged regions increases and 
expands, and on the other hand additional damaged areas arise, see 
Fig. 6.15c and Fig. 6.15d. At a certain loading level, a large share of 
the matrix is damaged, i.e., the density of microcracks is relatively 
high, compare Fig. 6.5d. Correspondingly, bundle damage evolves as 
the matrix damage increases. The sum of all these types of damage, 
matrix damage followed by bundle damage, eventually accumulates to 
macroscopic failure of the SMC composite. The full-field distribution of 
bundle damage is discussed in Appendix D. 


LR IN Dave à Q 1 
(a) Strain 1% (b) Strain 1.5 % (c) Strain 2% (d) Strain 2.5 % 


Figure 6.15: Evolution of matrix damage during increased loading in vertical direction, in 
analogy to experimental observations as presented in Fig. 6.5 


Our presented framework allows for the generation and computation 
of microstructures with orientations deviating from the planar isotropic 
state. The distribution of the predicted damage in the UPPH matrix due 
to a loading of 2.5 % strain in e,-direction is shown in Fig. 6.16. Basically, 
we observe a larger distribution of matrix damage for the planar isotropic 
orientation in comparison to orientations with a preferred alignment in 
loading direction. The less bundles are oriented in loading direction, 
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the less these bundles can bear the applied loading. Consequently, 
the loading level in the matrix is higher and so is the general damage 
level. Occurring damage is distributed in the transition areas between 
bundles at the bundle tips, and localizes especially at rims of matrix rich 
areas, see Fig. 6.16b and Fig. 6.16c. Corresponding analyzes regarding 
the distributions of bundle damage w.r.t. the orientation are given in 


Appendix D. 


(a) Planar isotropic orientation (b) Higher alignment in (c) Unidirectional orientation 
e,-direction in e„-direction 


Figure 6.16: Predicted UPPH matrix damage according to Sec. 6.3.3 for different 
orientations and loading of 2.5 % strain in ex-direction (see Fig. 6.12 for corresponding 
microstructures). Bundles are hidden. 


6.6.2 Structural analysis perspective 


Utilizing the experimental observations discussed in Sec. 6.2, we derive 
different macroscopic failure criteria based on corresponding full-field 
damage distributions. To investigate the influence of the interaction 
between fiber bundle orientation and loading direction on the damage 
evolution and the consequent point of failure, we load generated 
microstructures with different orientation states (see Fig. 6.12) in differ- 
ent, non-uniform distributed directions n € 8? = {x € R? | ||a|| = 1}. 
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Therefore, we sample 1000 directions on a unit half sphere with 
positive e,-direction, see Fig. 6.17. In analogy to the conducted 
experiments, we apply an effective uniaxial strain boundary condition 
in all sampled directions. 

For each direction, we compute the 
damage evolution and the stress 
response, and evaluate the macro- 
scopic stiffness reduction. For the 
latter, we compare the current effec- 
tive stiffness with the initial stiffness 
in terms of the directional Young’s 
modulus E(n) computed via the 


stress-strain relation Ao(n)/Ae(n). 


Given a direction n, the associated 
Figure 6.17: Sampled directions on unit 


directional stress a(n) can be ex- 
(a) half sphere 


tracted via 
o(n)=a-(n@n), (6.63) 


and similarly for the directional strain. 

Fig. 6.18 shows the resulting failure surfaces of a planar isotropic 
microstructure for different stiffness reductions, derived from the range 
of reductions known from the experimental observations, see Fig. 6.4d. 
These plots define the maximum allowable stresses for a given level of 
admissible damage before total failure. Hence, these failure surfaces 
serve as criterion on a macroscopic integration point level whether a 
part will sustain a certain stress level or not. 

For comparability, we depict the failure surfaces in terms of a constrained 
plot in the {e,,e,,e,}-space. Note that the negative axes for the stresses 
in the e,-e,-plane are caused by the corresponding loading directions 
that point into negative e,- and/or e,-direction, see Fig. 6.17. The 
associated directional stress levels themselves are always positive. 
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(a) 10 % reduction (b) 20 % reduction (c) 30 % reduction (d) 40 % reduction 


Figure 6.18: Failure surface plots for a microstructure with planar isotropic orientation 


Evaluating Fig. 6.18, the planar isotropic nature of the microstructure 
is apparent, as the failure surface is rotationally symmetric about the 
e,-axis, i. e., the directional stress evoking a certain stiffness reduction 
remains the same irrespective of the loading direction within the e,- 
e,-plane. Choosing an admissible stiffness reduction of 20% would 
allow for a directional stress of about 120 MPa in the e,-e,-plane, see 
Fig. 6.18b. An increasing e,-component of the loading direction leads to 
lower allowed stress levels before failure, as the supporting influence 
of the fiber bundles decreases for an increasing deviation of the loading 
direction away from the isotropic e,-e,-plane. 

Considering a microstructure with bundles being preferably aligned 
in e,-direction, see Fig. 6.12b, we obtain the failure surfaces shown in 
Fig. 6.19. The influence of the non-isotropic bundle orientation is notice- 
able through the elongated shape of the failure surface in e,-direction. 
As more bundles are aligned in e,-direction, the corresponding stress 
level is highest and decreases for loadings with increasing components 
in e,-direction and/or e,-direction. 

For an allowed stiffness reduction of 20 %, an SMC composite part with 
said orientation could bear stresses of up to 250 MPa in e,-direction and 
about 75 MPa in e,-direction, see Fig. 6.19b. If we also add a component 
in e,-direction with an angle of about 45°, the structure withstands 
stresses between 75 MPa and 105 MPa before total failure. 
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4007s 


(a) 10 % reduction (b) 20 % reduction (c) 30 % reduction (d) 40 % reduction 


Figure 6.19: Failure surface plots for a microstructure with higher alignment in 
e,-direction 


For a microstructure with (fully) aligned bundles in e,-direction, see 
Fig. 6.12c, the resulting failure surfaces follow a similar pattern as for 
the slight orientation in e,-direction shown in Fig. 6.19, but are more 
pronounced in their extremes, see Fig. 6.20. Furthermore, the failure 
surfaces are rotationally symmetric about the e,-axis, i. e., the admissible 
stress level remains unchanged for any loading direction in a e,-e,-plane. 
For all loading directions with no component in e.-direction, the allowed 
stress level is lowest. 


w 


(a) 10 % reduction (b) 20 % reduction (c) 30 % reduction (d) 40 % reduction 


Figure 6.20: Failure surface plots for a microstructure with unidirectional orientation in 
e,-direction 


The bulges, that become more pronounced for higher stiffness reductions, 
are caused by the reinforcing character of the fiber bundles for any 
loading direction with at least a small component in e,-direction. The 
bearable stress significantly increases for any loading direction with a 
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non-zero component in e,,-direction, see Fig. 6.20d. A further increase 
of the e,-component further increases the bearable stress level, but its 
effect decreases. 


6.6.3 Design perspective 


In addition to the structural analysis of our SMC composite, we use 
our approach to analyze the SMC composite loading conditions from 
a design perspective. For that case, we consider loadings that evoke 
specific stress levels, e. g., motivated by given component specifications 
or boundary conditions, and evaluate the resulting stiffness reduction. 
The gained knowledge on stiffness reduction due to given stress levels 
helps identifying and revisiting critical component areas during the 
design process to counter failure from the start. 

To scrutinize the relation between applied stress level and resulting 
stiffness reduction, we compute the directional-dependent relative 
residual Young’s modulus (RRYM) as 


_ Ein) 
Eo(n) 


ö(n) (6.64) 
via the current (damaged) Young’s modulus E(n) and the initial Young’s 
modulus Eo(n) in a direction n. Hence, a value of ö(n) = 100% 
represents a sound Young’s modulus in that direction n. 

Similar to the analysis conducted in Sec. 6.6.2, we depict the RRYM via 
iso-surface plots in the {e,, ey, e.}-space. Note that the negative axes for 
the RRYM in the e,-e,-plane are caused by the corresponding loading 
directions that point into negative e,- and/or e,-direction, see Fig. 6.17. 
The associated directional RRYM themselves are always positive. 

For a microstructure with planar isotropic orientation, the RRYM plots 
are shown in Fig. 6.21. For an applied stress of 10 MPa, the behavior 
is purely elastic and damage does not initiate yet, which yields a 
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value of ö(n) = 100% in all directions, see Fig. 6.21a. If the stress 
level is 60 MPa, damage evolves and the stiffness is reduced, see 
Fig. 6.21b. The stiffness reduction is smallest in the e,-e,-plane and 
highest for loading in out-of-plane e,-direction, which is in line with 
the observations made in Fig. 6.18. The RRYM is rotationally symmetric 
about the e,-axis due to the planar isotropic nature of the microstructure. 
For higher loadings, the RRYM decreases further, as shown in Fig. 6.21c. 
For convenience, we highlighted the undamaged regions with values 
of ö(n) = 100% in the e,-e,-plane in gray. A stress level of 180 MPa 
results in a RRYM within the e,-e,-plane of about 6(n) = 67 %, which 
corresponds to the experimental observations, compare Fig. 6.21d and 
Fig. 6.4. If we load the microstructure with a component in e,-direction 
having a 45° angle, the RRYM would be 6(n) = 53% for the same 
loading level. 
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(a) 10 MPa stress (b) 60 MPa stress (c) 120 MPa stress (d) 180 MPa stress 


Figure 6.21: Relative residual Young’s modulus plots for a microstructure with planar 
isotropic orientation 


The RRYM plots for an SMC composite microstructure that is preferably 
oriented in e,-direction are shown in Fig. 6.22. For stress levels 
that induce damage, the RRYM is highest in e,-direction due to the 
reinforcing nature of the bundles and decreases towards a loading in e,- 
direction. The RRYM is lowest for loadings in out-of-plane e.-direction, 
as we lack any supportive effect of the fiber bundles, see Fig. 6.22c. For 
an applied stress of 120 MPa, the RRYM is about 90 % in e,-direction and 
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about 60 % in e,-direction. Similar to observations made in Fig. 6.21d, 
we see the formation of bulges in Fig. 6.22d. Even a small proportion 
of bundles oriented in loading direction have a comparatively high 
influence on the overall damage evolution, and the RRYM in vicinity of 
the e,-direction immediately increases with a comparatively high slope. 
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(a) 10 MPa stress (b) 60 MPa stress (c) 120 MPa stress (d) 180 MPa stress 


Figure 6.22: Relative residual Young’s modulus plots for a microstructure with higher 
alignment in e,-direction 


The resulting RRYM plots for a microstructure with unidirectionally 
aligned bundles in e,.-direction are shown in Fig. 6.23. These plots are 
rotationally symmetric about the unidirectional e,-axis. The general 
behavior is similar to the one observed in Fig. 6.22, but more extreme. 
Again, we observe a higher sensitivity of the RRYM in vicinity to 
loadings in the e,-e,-plane and a lower sensitivity for loadings oriented 
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(a) 10 MPa stress (b) 60 MPa stress (c) 120 MPa stress (d) 180 MPa stress 


Figure 6.23: Relative residual Young’s modulus plots for a microstructure with unidirec- 
tional orientation in e,,-direction 
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6.7 Conclusions 


This chapter was devoted to identifying macroscopic and anisotropic 
failure criteria for SMC composites that are motivated by damage 
evolution on the microscale. Accumulating damage in matrix and 
bundles in terms of microcracking (Fitoussi et al., 2005; Ogihara 
and Koyanagi, 2010; Schöttl et al., 2021a) on the microscale yields a 
hardening regime on the macroscale and eventually ends in abrupt 
failure (Meraghni and Benzeggagh, 1995; Trauth, 2020). 


To capture such a structural behavior, we applied an anisotropic and 
modular framework (Görthofer et al., 2022b). The modular concept 
allows for the description of fully anisotropic stiffness degradation based 
on extraction tensors. Inspired by Puck’s theory for laminates (Puck 
and Schiirmann, 2002; Knops, 2008), as well as dilatation and distortion, 
we introduced a number of extraction tensors specifically designed to 
model damage evolution in SMC composites. 

With the help of Bayesian optimization (Frazier, 2018), integrating 
a Matern 5/2 kernel, an anisotropic distance function for the under- 
lying Gaussian process and an "expected improvement" acquisition 
function, we were able to identify all necessary parameters in the 
highly heterogeneous solution space. An implementation combining the 
experimental results, an FFT-based full-field homogenization approach 
(Schneider, 2021), generated SMC composite microstructures (Görthofer 
et al., 2020) and a corresponding error measure, allowed for a successful 
minimization of the cost function. 

We compared our predicted full-field damage results to corresponding 
in-situ CT scan analyzes. Matrix damage in the form of microcracks 
initiates at the rims of matrix rich regions and is followed by bundle 
damage. Utilizing generated SMC composite microstructures (Görthofer 
et al., 2020), we transfered our findings onto microstructures with 
individually selected orientations. 
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To propel scale-transition, we identified appropriate micromechanics- 
based anisotropic failure surfaces. We proposed, on the one hand, 
failure surfaces in stress space considering given admissible stiffness 
reductions, useful for structural analysis processes. On the other 
hand, we introduced relative residual stiffness (failure) surfaces based 
on predetermined component boundary conditions, which are vital 
for (virtual) design processes. We analyzed the influence of the 
fiber bundle orientation on both approaches for specific representative 
microstructures. 

The presented framework forms a foundation for a micromechanically 
and physically motivated database for anisotropic failure criteria on 
integration point level for SMC composites. In analogy to databases 
identified for plasticity and short fiber reinforced composites (Köbler 
et al., 2018; Gajek et al., 2021), a proper discretization and sampling of 
the orientation space can be established for orientation-dependent SMC 
composite failure criteria. Furthermore, an uncertainty quantification 
of the resulting damage and consequent macroscopic failure w.r.t. 
microscopic parameter variations (such as phase properties or volume 
fractions) can be investigated. The implemented damage model directly 
operates on the compliance tensor, which allows for a straightforward 
coupling of damage evolution to other phenomena such as plasticity 
(Simo and Ju, 1987; Ju, 1989) in future applications. 
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Chapter 7 


Summary and outlook 


In this thesis, we were concerned with full-field modeling of SMC 
composites. To provide an understanding of the SMC composite under 
consideration, we presented an overview on the manufacturing process 
and the resulting characteristic microstructure in Chapter 3. 

The main focus of our research regarded damage modeling of SMC 
composites on the microscale and a consecutive holistic approach to 
transfer the results to be used within macroscopic structural material 
models. The research aimed at a deeper understanding of the material 
behavior on the microscale to propel the use of SMC composites as 
structural, load-bearing components. Within this thesis, we focused on 
three main research topics. 

Having accurate SMC composite microstructures (including important 
properties such as orientation distribution) is a key foundation for 
corresponding material models to be developed, verified and validated. 
Hence, in Chapter 4, we presented a method to generate SMC composite 
microstructures using an adapted random sequential addition approach. 
As we relied upon an exact closure approximation in two dimensions 
and a quasi-random orientation sampling, we were able to generate 
high fidelity representative unit cells. Using our improved unit cell 
generator, we conducted a detailed sensitivity analysis via thousands 
of microstructures, i.e., we thoroughly investigated the influence of 
microstructural and mechanical parameters on the effective elastic 
properties of our SMC composite. These results, on the one hand 
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provide detailed relations between microscopic and elastic macroscopic 
parameters, and on the other hand serve as starting point for studying 
the inelastic behavior of SMC composites. 


In Chapter 5, we presented a general convex framework to compute 
anisotropic damage evolution via a direct operation on the compliance 
tensor. The model is formulated in the setting of generalized standard 
materials, is thermodynamically consistent and fulfills Wulfinghoff’s 
damage growth criterion. Using a modular ansatz based on extraction 
tensors, we were able to formulate the model so that it can be adapted 
and applied to any hardening-type material. Furthermore, we designed 
specific extraction tensors motivated by Puck’s theory on laminates in 
an averaged stress setting. A multitude of different examples showed 
the capability of predicting any anisotropic damage evolution on the 
micro- and the macroscale. In future applications, the modular character 
allows for an extension of the model to also account for softening-type 
behavior, as well as coupling to further phenomena such as plasticity. 

We applied our damage model to the considered SMC composite in 
Chapter 6. Upon analyzing the behavior of the material via different 
experimental investigations and in-situ CT testing in combination 
with results taken from corresponding publications, we identified 
the main damage phenomena on the microscale and the associated 
macroscopic structural behavior. We introduced new sets of extraction 
tensors accounting for damage in the fiber bundles and the matrix. 
Therefore, we used Puck-type criteria in combination with a pencil 
glide ansatz to account for maximum stresses, as well as ansatzes taking 
into account distortion and dilatation. A Bayesian optimization based 
on Gaussian distributions and an anisotropic distance function enabled 
identification of all associated damage parameters via corresponding 
comparisons to experimental results. Computations and FFT-based full- 
field homogenization on different generated microstructures loaded in 
1000 directions allowed for an holistic approach identifying macroscopic 
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failure criteria based on damage evolution on the microscale. Here, we 
provided both a design perspective and a structural analysis perspective 
in terms of different failure surfaces. A sampling of the orientation 
space for the presented approach can provide a database for orientation- 
dependent SMC composite failure criteria in future work and give a 
meaningful basis for macroscopic phenomenological failure modeling. 
Furthermore, an uncertainty quantification of microscopic anisotropic 
damage evolution and predicted macroscopic failure regarding varia- 
tions of microscopic properties can be conducted. 


Taking all results into consideration, this thesis provides tools to generate 
bundle microstructures, model any hardening-type damage material 
and a corresponding setup to identify all necessary damage parameters 
in the framework of FFT-based full-field computations. Furthermore, 
we applied our models to account for damage evolution in SMC 
composites on the microscale in combination with an approach to 
identify anisotropic failure criteria on the macroscale. Our research can 
be extended to describe SMC composites in more detail via an adaption 
of the bundle generator or the damage model via further experimental 
investigations (Schemmann et al., 2018c;a). Furthermore, an extension 
towards fatigue (Bartkowiak et al., 2020; Magino et al., 2022) would be 
possible. A coupling to sophisticated compression molding flow models 
(Meyer et al., 2020) could form the basis of an improved virtual process 
chain for SMC composites (Görthofer et al., 2019b). Via a combination 
with a proper plasticity model and a new set of extraction tensors, the 
damage model could be used to describe DiCoFRTP such as LFT in the 
third doctoral generation of IRTG 2078. 
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Appendix A 


Proofs for the explicit exact 
closure approximation 


First, we wish to show the first identity of (3.13). (The second identity 
follows by symmetry considerations.) Thus, we wish to prove the 


equation 
2T 2 
1 cos“ & ge 1 l (A.1) 
2n Jo By cos? d+ Basin“ d 1+ By 
Denote the left hand side of (A.1) by Aı. By symmetry, we have 
2 3 1 By 
A = d th a= =. A.2 
A =) 1+ atan? ¢ p A Bı (A.2) 


Substituting x = tan ¢ yields 


2 = 1 
= dz. A. 
Ai n Bı | (1+ aa?)(1 + z?) a ee) 


Decomposing the integrand into partial fractions yields’ 


1 ad 1 1 
(l1+aa?)1+2?) a-Iil+ar? a-—1142? 


(A.4) 


1 Ifa = 1, temporarily set a = 1 — e for e > 0 small, and let e — 0 in the end. 
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Thus, 


2 a = 1 1 == 
A = d dz}. (A5 
1 By lf Ita” —/ 1+? | (aa 


Using the identities (easily shown by substitution) 


= 1 1 % l 

= d A.6 
/ Eu. — | ite’ ee 
wel T 

Der A.7 
fae n (A7) 


we arrive at the identity 


= lya-1i 1 1 1 1 DR 
 Bıa-1 Biva+1 YByB+vB 1+Bi 


where we have used B, Bz = 1 for the last equality. Thus, we have 
shown (A.1). 


Next, we wish to show Eq. (3.14)ı. By definition we have 


i ge cos* & 
A = dd. AI 
ern = Bı cos? 6 + Basin? & o ee) 


Ay 


(A.8) 


Following identical steps as for A; we arrive at the identity 


2 = 1 
A = x. Al 
By f (1 + ax?)(1 + 22)? ee N) 


Decomposing the integrand into partial fractions yields 
1 — 
(1+ ax?)(1+22)2 


a? 1 a 1 is 1 1 
(l-a)? 1+ar? (l-a)? 1+r? 1-a (1+ zr?) 


(A.11) 
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Notice 
f > 1 — dx = a (A.12) 
3 — dx = F (A.13) 
m a oe de =. (A.14) 


The a u can be evaluated, for instance, by differentiating 
I =. de = I Z at z = 1. Consequently, we get 


_ 1 fa(va-1) 1 1 
Ayu = B; | (l-g ! | (A.15) 


which we rewrite, using the second identity in (A.8), 


FRE 1 1 2a(ya-l) | 
HH 2B (1+ va) (1 -= Va) lza ` 
(A.16) 
_ Aı 1 2a 
2 |1-ya 1-e|' 
The last factor can be further processed, 
1 2a 1+ vya -— 2a a— ya Ja 
1-ya 1-a l-a l-a 1+ vya 
VE; ‘ (A.17) 
= 2 _=1+ =1+A,, 


"Br Bi Bı +1 


where for the last two transformations we used Bı B2 = 1 and (A.8). 
Thus, we finally arrive at 


1 
Aun = z4 + Ai), (A.18) 


which we wanted to prove. 
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Equation (3.14) follows similarly to Eq. (3.14);. Equation (3.14)4 and 
Eq. (3.14); are trivially satisfied as the integrand is an odd function 
(w.r.t. m). Last but not least, Eq. (3.14)3 follows from the identity 


Ay = Ajit + 41122 (A.19) 


which can be either proven abstractly, i. e., Ajj = Mii Aijkk holds, or 
concretely for the integrals using cos? ¢ + sin? ¢ = 1. 
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Appendix B 


Evolution equations for internal 
variables 


In this appendix, we give a short overview on the formulations 
regarding the evolution of the internal variables and the associated 
Karush-Kuhn-Tucker (KKT) conditions of our damage model as 
presented in Sec. 5.2. 


Summary of formulations for evolution equations for 


internal variables of the damage model presented in Sec. 5.2 


General formulation z=(S,q) 659) 


D*+ hi 6.20) 
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Summary of Karush-Kuhn-Tucker (KKT) conditions for 


internal variables of the damage model presented in Sec. 5.2 


Karush-Kuhn-Tucker (KKT) conditions (5.23) 


420, file,a) <0, di filo, qi) = 
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Appendix C 


Convexity of the Hookean elastic 
energy density 


In this appendix, we compute the Hessian of the Hookean elastic energy 
(density) we, see formula (5.11). We show the Hookean elastic energy 
(density) to be jointly convex in both the strain and the compliance. This 
characteristic feature in combination with a hardening-type damage 
potential enables us to design a damage evolution framework which is 
convex itself. 

For a twice (Fréchet) differentiable function f : U C X — Ronan open 
subset of a (Banach) vector space, the Hessian at some point x € U may 
be represented as a quadratic form 


D? f(z): X >R, (C.1) 
which may be computed by the directional second derivative 


1 qd? 


D?’ f(x)[y] = apo + ày) N (C.2) 


For our problem at hand, we have U = Sym(d)x Sq and 

X = Sym(d) x Sym(Sym(d)). Then, the Hessian of the elastic energy 
(5.11) 

at S dx S R 

w ym(d) x Sa > 7 (C3) 

we(e,S) = ze-S”|e 
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at (e, S) in direction (€,L) computes as 


1 d? /1 
Duele SEL gr (3E4409: 6407 etaa) | 
1 d 


A=0 


1 —1 =] 1 -4 
ze SLS [e] + 36: S7 [E] 
ze STILS? [6] + le “SLS LS" [e] 
E- S7? [E] — €- S LS- [e] + Le -S“LS“ LS“? [e] 


(€ - LS“! [e]) -S~! [€ - LS“ [el]. 
(C.4) 


We see that for any admissible strain and compliance 
(e,S) € Sym(d) x Sg and any direction (&,L) € Sym(d) x Sym(Sym(d)), 
the Hessian is non-negative and hence the elastic energy is convex. 
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Appendix D 


Full-field damage evolution for 
different microstructures 


An increase of the total loading in combination with evolving matrix 
damage leads to an increase of load being distributed onto specific fiber 
bundles. Inevitable, damage is also evoked within these bundles in the 
form of microcracks running through the bundles along the principal 
direction and interface debonding. 

The evolution of bundle damage associated to normal stresses perpen- 
dicular to the bundle direction is shown in Fig. D.1 for different loading 
levels in vertical direction (in analogy to Fig. 6.15 and Fig. 6.5). The 
higher the loading, the higher the number of bundles that are damaged 
and the higher the damage level in certain bundles. Generally, bundles 
perpendicular to the loading direction encounter the highest damage, 
whereas bundles in loading direction are hardly damaged for the case at 
hand, see Fig. D.1d. 

Analogously, we encounter damage in bundles due to shear in lon- 
gitudinal bundle direction, which is shown in Fig. D.2. For higher 
loading levels, bundles that a preferably oriented in loading direction are 
damaged according to the Puck-type shear criterion, whereas bundles 
perpendicular to the loading direction are hardly damaged. Hence, 
our applied Puck-type cases as discussed in Sec. 6.3.4 and Sec. 6.3.4 
are somewhat complementary in their effect on the overall damage 
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Figure D.1: Evolution of bundle damage perpendicular to bundle direction during 
increased loading in vertical direction, in analogy to experimental observation as presented 
in Fig. 6.5 


evolution, offering a full coverage of possible stiffness degradations in 
the bundles. 
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Figure D.2: Evolution of bundle damage in bundle direction during increased loading in 
vertical direction, in analogy to experimental observation as presented in Fig. 6.5 


The evolution and distribution of matrix damage is affected by the 
orientation of the SMC composite microstructure, see Sec. 6.6.1. The 
alignment of bundles is essential, as their reinforcing character is more 
pronounced the more bundles point into loading direction. Corre- 
spondingly, damage in bundles due to normal stresses perpendicular 
to the principal directions is higher and further distributed, the less 
fibers are oriented in loading direction, see Fig. D.3. In the limit of 
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a pure unidirectional orientation, bundles do not undergo damage 


perpendicular to their principal direction, see Fig. D.3c. 
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Figure D.3: Predicted bundle damage according to Sec. 6.3.4 for different orientations and 
loading in e,-direction (see Fig. 6.12 for corresponding microstructures). Matrix is hidden. 


Damage due to shear stresses in principal bundle direction is also 
affected by the orientation of the SMC composite microstructure, see 
Fig. D.4. In concurrence with the other damage cases (see Fig. 6.16 and 
Fig. D.3), the damage level of bundles for the considered case is higher, 
if more bundles are aligned in loading direction, see Fig. D.4b. For an 
orientation state in vicinity to a unidirectional alignment in loading 
direction, damage due to shear stresses localizes mainly at bundle tips 
at the rims of matrix rich areas, see Fig. D.4c. In these regions, bundles 
are subjected to local stress excesses due to the prior matrix damage, see 
Fig. 6.16c, and hence damage evolution is accelerated. 
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(a) Planar isotropic orientation (b) Higher alignment in es- (c) Unidirectional orientation 
direction in e,-direction 


Figure D.4: Predicted bundle damage according to Sec. 6.3.4 for different orientations and 
loading in e-direction (see Fig. 6.12 for corresponding microstructures). Matrix is hidden. 
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